# A generalized curvilinear solver for spherical shell Rayleigh-Bénard
convection

Souvik Naskar  Karu Chongsiripinyo  Anikesh Pal<EMAIL_ADDRESS>Akshay
Jananan  Department of Mechanical Engineering, Indian Institute of Technology,
Kanpur 208016, India  Department of Mechanical Engineering, Chulalongkorn
University, Bangkok, Thailand

###### Abstract

A three-dimensional finite-difference solver has been developed and
implemented for Boussinesq convection in a spherical shell. The solver
transforms any complex curvilinear domain into an equivalent Cartesian domain
using Jacobi transformation and solves the governing equations in the latter.
This feature enables the solver to account for the effects of the non-
spherical shape of the convective regions of planets and stars. Apart from
parallelization using MPI, implicit treatment of the viscous terms using a
pipeline alternating direction implicit scheme and HYPRE multigrid accelerator
for pressure correction makes the solver efficient for high-fidelity direct
numerical simulations. We have performed simulations of Rayleigh-Bénard
convection at three Rayleigh numbers $Ra=10^{5},10^{7}$ and $10^{8}$ while
keeping the Prandtl number fixed at unity ($Pr=1$). The average radial
temperature profile and the Nusselt number match very well, both qualitatively
and quantitatively, with the existing literature. Closure of the turbulent
kinetic energy budget, apart from the relative magnitude of the grid spacing
compared to the local Kolmogorov scales, assures sufficient spatial
resolution.

††journal: Journal of Computational Physics

##  1 Introduction

Fig. 1: Spherical shell geometry.

Turbulent thermal convection is ubiquitous in nature as a primary driving
mechanism for atmospheric and oceanic circulations [1]. Such convective
motions in Earth’s outer core or in the solar convective zone, for example,
provide energy to sustain global-scale magnetic fields in planets and stars[2,
3]. Such flow phenomena are further enriched due to the presence of global
rotation, external or self-generated magnetic fields, chemical reactions,
phase change, the porosity of the medium, and particle suspension [4].
Furthermore, the design of heat exchangers, cooling systems for electronics,
and indoor air circulation systems requires a fundamental understanding of
thermal convection [5, 6]. Rayleigh-Bénard convection (RBC) is a simple model
of thermal convection, where a fluid layer between two parallel plates is
heated from below and cooled from above. Such a plane layer geometry can be
considered, for example, as a local approximation of the tangent cylinder
region of Earth’s outer core, which is situated between the top and bottom
surfaces of the solid inner core and extending towards the north and south
poles, respectively, up to the core-mantle boundary. Simulations in the plane
layer geometry can reproduce the basic force balance and heat transfer
behavior that can be validated from well-designed laboratory experiments.
Therefore, this flow configuration has been extensively studied, with the
individual or combined effect of global rotation and magnetic fields [7, 8] to
model various geophysical and astrophysical turbulent flows [9].

In the geophysical and astrophysical context, however, a spherical shell
geometry is more pertinent to modeling planetary cores or stellar convective
zones. The most extensive body of literature in this geometry focuses on
”geodynamo” simulations that attempt to model convection in Earth’s outer core
convection and the associated geomagnetic field originating from it [10].
Mantle convection [11], rapidly rotating convection [12, 13], RBC without
rotation and magnetic field [14, 15, 16], deep convection in gas giants[17,
18], and solar convection [19] are among the other prolific areas of research
where spherical shell models are implemented. The superiority of these models
lies in their capability to model many essential dynamical features of
planetary atmospheres, such as thermal winds, strong shear layers, magnetic
buoyancy, meridional circulations, and large-scale flows. They can also
incorporate important geometric constraints, such as tangent cylinders and
curvature effects near the boundaries, whose combined or individual influence
can not be accounted for in a local Cartesian plane layer configuration [20].

The local plane layer and the global spherical shell simulations differ
primarily in the direction of gravity, which is generally kept vertically
downwards in the local Cartesian models. In contrast, the direction is
radially inwards in global spherical shell models. Additionally, rotating
convection in spherical shells exhibits distinct scales in the radial, axial,
and azimuthal directions [21], whereas, for the local Cartesian model, we need
to consider only two spatial scales: the horizontal scale of convection and
the vertical scale over which convection occurs. For both geometries, the
governing non-dimensional parameters are the Rayleigh numbers ($Ra$), which is
a non-dimensional measure of the thermal forcing, and the Prandtl number
($Pr$), representing the viscous to thermal diffusivity ratio. Apart from
this, the flow properties may also depend on the aspect ratio $\Gamma=W/H$
(where W and H are the horizontal and vertical extents of the domain) and the
radius ratio $\Gamma=r_{i}/r_{o}$ in-plane layer and spherical shell
geometries, respectively. The important global diagnostic quantities are the
Nusselt number $Nu$ and the Reynolds number $Re$, representing the non-
dimensional heat transfer and flow speed.

An intriguing question in this research direction is the scaling relation
between such a diagnostic quantity with a governing input parameter, such as
$Nu$, as a function of $Ra$. The thermal convection in planets and stars
occurs at parameter values that are several orders of magnitude away from the
reach of state-of-the-art numerical simulations and experiments. Therefore,
these scaling relations are valuable tools to extrapolate the results of these
experiments and simulations to planetary and stellar convective regimes. For
plane layer geometry, a Nusselt number scaling of $Nu\thicksim Ra^{2/7}$ is
found for moderate thermal forcing ($Ra\leqslant 10^{10}$), whereas, for
higher thermal forcing, a scaling relation of $Nu\thicksim Ra^{1/3}$ has been
widely reported [22]. A systematic investigation has been reported by [14],
who found the same scaling laws for the Nusselt number in the spherical
geometry. It should be noted here that though the global diagnostic quantities
exhibit similar behaviour, the local properties, such as the thickness of the
viscous and thermal boundary layers, are markedly different in the two
geometries. For example, the effect of curvature and a radially varying
gravitational acceleration (as appropriate in Earth’s core) results in
asymmetric boundary layers in the spherical geometry, in contrast to the
symmetric boundary layers in a plane layer geometry.

Experimental difficulties related to the radial direction of gravity make the
advances in spherical shell convection almost entirely dependent on massively
parallel numerical simulations. Existing solvers [23],[15] use spherical
harmonic decomposition of the flow variables in the azimuthal and latitudinal
directions while Chebyshev polynomials are used in the radial direction for
proper resolution of the boundary layers. In this paper, we report on the
development, implementation, and validation of a new finite-difference solver
for studying spherical shell convection. The solver can map any three-
dimensional curvilinear geometry to a computational Cartesian domain using the
Jacobi transformation. This enables us to solve the conservation equations in
Cartesian coordinates, which are much simpler than their spherical coordinate
counterpart, even after their modification by the Jacobi, elongation, and
stiffness matrix coefficients. Furthermore, the effect of the ellipticity of
the core-mantle boundary [24] and the anisotropic shape of the inner core [25]
on the azimuthal and latitudinal variation of radial heat flux can be
accounted for. The capability to account for any effect of the non-spherical
boundaries is the primary motivation for developing the present solver. The
solver uses second-order central spatial discretization, while temporal
discretization is achieved with the fractional step method [26]. In order to
avoid the stiffness induced by the fine resolution near the boundary layers,
the viscous terms have been treated implicitly, while the other terms are
marched explicitly. The fractional step marches the velocity field into an
intermediate field by a combination of the Alternating Direction Implicit
method (ADI), the Crank-Nicolson method (CN), and the third-order low-storage
Runge-Kutta method (RKW3) [26]. The remaining procedure in the fractional step
method is to remove the divergence residual from the velocity field after the
end of each RKW3 step, which in turn is achieved by pressure correction. We
use the multigrid HYPRE module to accelerate the pressure correction. The rest
of the article is structured as follows. Section 2 discusses the governing
equation used. The numerical scheme is described in 3. Results are presented
in Section 4 and summarized in section 5.

##  2 Governing Equations

We aim to investigate Rayleigh-Bénard convection of an incompressible,
Newtonian, Boussinesq fluid in a spherical shell geometry as illustrated in
figure 1. The spherical shell has an inner radius $r_{i}$ and an outer radius
$r_{o}$ kept at constant temperatures $T_{i}$ and $T_{o}$, respectively. The
shell gap $d=r_{o}-r_{i}$, the temperature difference $\Delta T=T_{i}-T_{o}$,
and the free-fall velocity $u_{f}=\\{g_{0}\alpha(T_{i}-T_{o})d\\}^{1/2}$ have
been used as the characteristics scale for length, temperature, and velocity,
respectively, to nondimensionalize the governing equations. Here, $g_{0}$ is
the gravitational acceleration at the outer radius. The relevant fluid
properties are the kinematic viscosity ($\nu$), thermal diffusivity
($\kappa$), and thermal expansion coefficient ($\alpha$). The non-dimensional
governing equations are expressed below using a Cartesian coordinate system.

$\frac{\partial u_{j}}{\partial x_{j}}=0,$ (1) $\frac{\partial u_{i}}{\partial
t}+\frac{\partial u_{i}u_{j}}{\partial x_{j}}=-\frac{\partial p}{\partial
x_{i}}+gT\delta_{in}+\sqrt{\frac{Pr}{Ra}}\frac{\partial^{2}u_{i}}{\partial
x_{j}\partial x_{j}},$ (2) $\frac{\partial T}{\partial t}+u_{j}\frac{\partial
T}{\partial x_{j}}=\frac{1}{\sqrt{RaPr}}\frac{\partial^{2}T}{\partial
x_{j}\partial x_{j}},$ (3)

where $g=(r_{o}/r)^{2}$ is the radial variation of gravitational acceleration
and
$\delta_{in}=\cos\theta\;\delta_{i1}+\sin\theta\cos\phi\;\delta_{i2}+\sin\theta\sin\phi\;\delta_{i3}$.
Here $\theta$ and $\phi$ are the colatitude and longitude as shown in figure
1. The non-dimensional temperature difference is defined as
$T=(T_{f}-T_{o})/(T_{i}-T_{o})$, where $T_{f}$ is the temperature of the
fluid. The non-dimensional parameters in these equations are the Rayleigh
number and the Prandtl number defined below.

$Ra=\frac{g_{0}\alpha\Delta Td^{3}}{\kappa\nu},Pr=\frac{\nu}{\kappa}$ (4)

In the subsequent section, we will use a coordinate transformation to convert
the spherical domain to a Cartesian domain.

##  3 Numerical Algorithms

###  3.1 Coordinate Transformation

(a) (b)

Fig. 2: (a) Physical curvilinear domain depicting the left half of the
$x_{3}=0$ plane in the spherical geometry with the corresponding (b)
transformed Cartesian computational domain.

To solve the governing equations 1-3 in a generalized curvilinear coordinate
system we perform coordinate transformation. The basic idea behind a
coordinate transformation is to transform a set of physical laws written in
Cartesian coordinates $x_{1},x_{2},x_{3}$ into an alternative form based on
generalized curvilinear coordinates $\zeta,\eta,\xi$ [26].

Physical law written in the Cartesian system Physical curvilinear grid

$\Bigg{\downarrow}\;\text{Grid transformation}$

Physical law written in the generalized system Computational Cartesian grid,
Jacobi terms

Such a transformation will result in the inclusion of additional coefficients
in the space derivatives in the governing equations, and the relation of this
transformation between the Cartesian and the generalized curvilinear
coordinate system is stored in a Jacobi matrix ($J$). The continuity,
momentum, and energy equations after the transformation are expressed below.

$\frac{\partial\left[C_{nj}u_{j}\right]}{\partial\zeta_{n}}=0$ (5)
$\frac{\partial\left|J^{-1}\right|u_{i}}{\partial
t}+\frac{\partial\left[C_{nj}u_{j}\right]u_{i}}{\partial\zeta_{n}}=-\frac{\partial
C_{ni}P}{\partial\zeta_{n}}+gT\delta_{in}+\sqrt{\frac{Pr}{Ra}}\frac{\partial}{\partial\zeta_{n}}\left(G_{nj}\frac{\partial
u_{i}}{\partial\zeta_{j}}\right)$ (6) $\frac{\partial T}{\partial
t}+\frac{\partial\left[C_{nj}u_{j}\right]T}{\partial\zeta_{n}}=\frac{1}{\sqrt{RaPr}}\frac{\partial}{\partial\zeta_{n}}\left(G_{nj}\frac{\partial
T}{\partial\zeta_{j}}\right)$ (7)

Here, $x_{i}$ denotes the coordinate $i$ of the Cartesian system, and
$\zeta_{i}$ denotes the coordinate $i$ of the generalized system. The
notations $x_{i}=(x_{1},x_{2},x_{3})=(x,y,z)$ and
$\zeta_{i}=(\zeta,\eta,\xi)=(\zeta_{1},\zeta_{2},\zeta_{3})$ have been used
interchangeably. After the transformation, the transformed governing equations
are solved as if in a Cartesian system. In this context, grid transformation
is often synonymously used with coordinate transformation as the curvilinear
domain (i.e., a spherical shell domain in our case) is transformed into a new
computational Cartesian domain. Here _$J^{-1}$ ,$C_{ij}$_ and _$G_{ij}$_ are

$J^{-1}=\left[\begin{array}[]{lll}\partial x_{1}/\partial\zeta&\partial
x_{2}/\partial\zeta&\partial x_{3}/\partial\zeta\\\ \partial
x_{1}/\partial\eta&\partial x_{2}/\partial\eta&\partial x_{3}/\partial\eta\\\
\partial x_{1}/\partial\xi&\partial x_{2}/\partial\xi&\partial
x_{3}/\partial\xi\end{array}\right]=:\left[\partial
x_{i}/\partial\zeta_{j}\right]$ (8)

$C_{ij}=\left|J^{-1}\right|\frac{\partial\zeta_{i}}{\partial x_{j}}\quad
G_{ij}=\left|J^{-1}\right|\frac{\partial\zeta_{i}}{\partial
x_{k}}\frac{\partial\zeta_{j}}{\partial x_{k}}$ (9)

The determinant _$|J^{-1}|$_ is the volume ratio of the original cell to the
transformed cell, whereas _$C_{ij}$_ and _$G_{ij}$_ are grid elongation and
skewness coefficients, respectively. The side length, and consequently the
side area and the total volume, of a transformed cell is chosen to be unity.

###  3.2 Jacobi terms

Fig. 3:  A transformed computational cell associated with the grid point
$(i,j,k)$; where $i$, $j$, $k$ are the integer indices used to identify
discrete space in the $\zeta$, $\eta$, and $\xi$ directions, respectively.

Figure 3 demonstrates a cell $(i,j,k)$ in a transformed computational domain.
The Jacobi terms, $J^{-1}$, $C_{pq}$, and $G_{pq}$, as expressed in equations
8 and 9 are stored at the cell’s faces. The calculation of $J^{-1}$ , $C_{pq}$
and $G_{pq}$ is given below.

  1. 1.

$J^{-1}$ is computed at every cell face, denoted by $J^{-1,fc}$; where $fc$
indicates cell face (1-3), by calculating all the nine components in $J^{-1}$.
For instance, we can compute the components of $J^{-1,2}$ of a cell (i,j,k) as
follows,

$\displaystyle(\partial\overrightarrow{x}/\partial\zeta)|_{i,j,k}=0.125*($
$\displaystyle+\overrightarrow{x}|_{i+1,j+1,k}+\overrightarrow{x}|_{i+1,j,k}+\overrightarrow{x}|_{i+1,j+1,k+1}+\overrightarrow{x}|_{i+1,j,k+1}$
$\displaystyle-\overrightarrow{x}|_{i-1,j+1,k}-\overrightarrow{x}|_{i-1,j,k}-\overrightarrow{x}|_{i-1,j+1,k+1}-\overrightarrow{x}|_{i-1,j,k+1})$
$\displaystyle(\partial\overrightarrow{x}/\partial\eta)|_{i,j,k}=0.5*($
$\displaystyle+\overrightarrow{x}|_{i,j+1,k}+\overrightarrow{x}|_{i,j+1,k+1}-\overrightarrow{x}|_{i,j,k}-\overrightarrow{x}|_{i,j,k+1})$
$\displaystyle(\partial\overrightarrow{x}/\partial\xi)|_{i,j,k}=0.5*($
$\displaystyle+\overrightarrow{x}|_{i,j+1,k+1}+\overrightarrow{x}|_{i,j,k+1}-\overrightarrow{x}|_{i,j,k}-\overrightarrow{x}|_{i,j+1,k}).$

  2. 2.

Calculate $det(J^{-1,fc})$, denoted by $|J^{-1,fc}|$.

  3. 3.

The variable $|J^{-1}|$ in equation 6 is an averaged value at the cell center
calculated from the six surrounding faces:

$|J^{-1}|_{i,j,k}=\frac{1}{6}\left(\sum_{fc=1}^{3}|J^{-1,fc}|_{i,j,k}+|J^{-1,1}|_{i,j+1,k}+|J^{-1,2}|_{i+1,j,k}+|J^{-1,3}|_{i,j,k+1}\right).$

  4. 4.

Compute $J^{fc}=\left[\partial\zeta_{i}/\partial x_{j}\right]$ simply by the
straight-forward inversion, $J^{-1,fc}$:
$\\{J^{-1,fc}\\}^{-1}=det(J^{-1,fc})^{-1}\\{cof(J^{-1})\\}^{T}$

  5. 5.

Calculate $C_{pq}$ and $G_{pq}$ at face $fc$, denoted by $C^{fc}_{pq}$ and
$G^{fc}_{pq}$ from $J^{fc}$ using equation 9.

###  3.3 Spatial discretization

Fig. 4:  Stencils used for computing 11, 12, and 13.

The spatial derivatives in 6 and 7 are discretized using a second-order
central finite difference scheme. Figure 4 illustrates the stencils used to
discretize the term 10 using this scheme.

$\frac{\partial}{\partial\zeta_{p}}\left[G_{pq}\frac{\partial\phi}{\partial\zeta_{q}}\right]$
(10)

Equation 10 consists of $9$ terms. We present the discretization of term 1
($p$=1 and $q$=1), term 2 ($p$=1 and $q$=2), and term 3 ($p$=1 and $q=3$) as
examples.

$\displaystyle\left(\frac{\delta}{\delta\zeta_{1}}\left[G_{11}\frac{\delta\phi}{\delta\zeta_{1}}\right]\right)_{i,j,k}=$
$\displaystyle\left[G_{11}\frac{\delta\phi}{\delta\zeta_{1}}\right]_{i+1/2,j,k}-\left[G_{11}\frac{\delta\phi}{\delta\zeta_{1}}\right]_{i-1/2,j,k}$
$\displaystyle=$
$\displaystyle+G^{2}_{11}|_{i+1,j,k}[\phi|_{i+1,j,k}-\phi|_{i,j,k}]$
$\displaystyle-G^{2}_{11}|_{i,j,k}\quad[\phi|_{i,j,k}-\phi|_{i-1,j,k}]$ (11)
$\displaystyle\left(\frac{\delta}{\delta\zeta_{1}}\left[G_{12}\frac{\delta\phi}{\delta\zeta_{2}}\right]\right)_{i,j,k}=$
$\displaystyle\left[G_{12}\frac{\delta\phi}{\delta\zeta_{2}}\right]_{i+1/2,j,k}-\left[G_{12}\frac{\delta\phi}{\delta\zeta_{2}}\right]_{i-1/2,j,k}$
$\displaystyle=$
$\displaystyle+\frac{G^{2}_{12}|_{i+1,j,k}}{2}\left[+\frac{\phi|_{i,j+1,k}+\phi|_{i+1,j+1,k}}{2}\right.$
$\displaystyle\qquad\qquad\qquad\left.-\frac{\phi|_{i,j-1,k}+\phi|_{i+1,j-1,k}}{2}\right]$
$\displaystyle-\frac{G^{2}_{12}|_{i,j,k}}{2}\quad\left[+\frac{\phi|_{i,j+1,k}+\phi|_{i-1,j+1,k}}{2}\right.$
$\displaystyle\quad\qquad\qquad\>\>\>\>\,\left.-\frac{\phi|_{i,j-1,k}+\phi|_{i-1,j-1,k}}{2}\right]$
(12)
$\displaystyle\left(\frac{\delta}{\delta\zeta_{1}}\left[G_{13}\frac{\delta\phi}{\delta\zeta_{3}}\right]\right)_{i,j,k}=$
$\displaystyle\left[G_{13}\frac{\delta\phi}{\delta\zeta_{3}}\right]_{i+1/2,j,k}-\left[G_{13}\frac{\delta\phi}{\delta\zeta_{3}}\right]_{i-1/2,j,k}$
$\displaystyle=$
$\displaystyle+\frac{G^{2}_{13}|_{i+1,j,k}}{2}\left[+\frac{\phi|_{i,j,k+1}+\phi|_{i+1,j,k+1}}{2}\right.$
$\displaystyle\qquad\qquad\qquad\left.-\frac{\phi|_{i,j,k-1}+\phi|_{i+1,j,k-1}}{2}\right]$
$\displaystyle-\frac{G^{2}_{13}|_{i,j,k}}{2}\quad\left[+\frac{\phi|_{i,j,k+1}+\phi|_{i-1,j,k+1}}{2}\right.$
$\displaystyle\quad\qquad\qquad\>\>\>\>\,\left.-\frac{\phi|_{i,j,k-1}+\phi|_{i-1,j,k-1}}{2}\right]$
(13)

###  3.4 Temporal discretization

For temporal discretization, a fractional step method is used where a velocity
field is sequentially advanced in multiple substeps. We use a combination of
the Alternating Direction Implicit method (ADI), the Crank-Nicolson method
(CN), and the third-order low-storage Runge-Kutta method (RKW3) to march to an
intermediate field as described below [26].

####  3.4.1 Alternating Direction Implicit method

Alternating Direction Implicit (ADI) method has been used to treat the viscous
term implicitly while marching in one direction at a time. We demonstrate the
method with a two-dimensional diffusion equation as shown in equation 14. To
solve 14 using the Euler method, the procedure is to perform implicit Euler in
the $x$ direction with explicit Euler in the $y$ direction for the first half
($\Delta t/2$), and vice versa for the second half $\Delta t/2$ as shown in 15
and 16.

$\frac{\partial\phi}{\partial t}=\alpha\left[\frac{\partial^{2}\phi}{\partial
x^{2}}+\frac{\partial^{2}\phi}{\partial y^{2}}\right]$ (14)

$\frac{\phi^{n+\frac{1}{2}}-\phi^{n}}{\Delta
t/2}=\alpha\left[\frac{\partial^{2}\phi^{n+\frac{1}{2}}}{\partial
x^{2}}+\frac{\partial^{2}\phi^{n}}{\partial y^{2}}\right]$ (15)

$\frac{\phi^{n+1}-\phi^{n+\frac{1}{2}}}{\Delta
t/2}=\alpha\left[\frac{\partial^{2}\phi^{n+\frac{1}{2}}}{\partial
x^{2}}+\frac{\partial^{2}\phi^{n+1}}{\partial y^{2}}\right]$ (16)

####  3.4.2 Crank-Nicolson method

The Crank-Nicolson (CN) method splits the right-hand side into two equal
parts, the implicit and the explicit, as demonstrated in 17 and 18.

$\frac{\partial\phi}{\partial t}=\alpha\frac{\partial^{2}\phi}{\partial
x^{2}}$ (17)

$\frac{\phi^{n+1}-\phi^{n}}{\Delta
t}=\frac{\alpha}{2}\left[\frac{\partial^{2}\phi^{n+1}}{\partial
x^{2}}+\frac{\partial^{2}\phi^{n}}{\partial x^{2}}\right]$ (18)

####  3.4.3 Third order Runge-Kutta method

Substep |  $h$ |  $\beta$ |  $\Pi$
---|---|---|---
1 | 8$\Delta t$/15 | 1 | 0
2 | 2$\Delta t$/15 | 25/8 | -17/8
3 | 1$\Delta t$/3 | 9/4 | -5/4
Table 1: RKW3 parameters.

The third-order low-storage Runge-Kutta method (RKW3) uses only two storage
variables. Marching is accomplished in three substeps, briefly summarized
here. Given an equation for $\phi$,

$\frac{\partial\phi}{\partial t}=\text{R}(\phi).$ (19)

RKW3 is implemented in the following manner,

$\frac{\phi^{rk}-\phi^{rk-1}}{{h^{rk}}}={\beta^{rk}}\text{R}(\phi^{rk-1})+{\Pi^{rk}}\text{R}(\phi^{rk-2})\,.$
(20)

Here, $rk$ goes from substep 1 to substep 3, and the values of ${h}$,
${\beta}$, and ${\Pi}$ are given in table 1.

####  3.4.4 The ADI-CN-RKW3 combined marching scheme

Fig. 5: ADI-CN-RKW3 combined marching scheme. PC denotes the pressure
correction.

The above-mentioned algorithms (ADI, CN, and RKW3) are combined to march the
governing equations to an intermediate state temporally. The right-hand side
of equation 6 is split into explicit and implicit terms as indicated by the
subscripts $ex$ and $im$ in equation 21. Depending on the grid skewness
$G_{ij}$, the diagonal components of the viscous terms are susceptible to the
stiffness of the discretized systems and are, therefore, marched implicitly.
The ADI scheme is used since there are three viscous terms containing
$G_{11}$, $G_{22}$, and $G_{33}$. At a given time, they are split into two
parts using CN. These steps are shown in 22, 23, and 24 as an example for
substep 1 of the RKW3 marching scheme.

$\displaystyle\begin{split}\frac{\partial J^{-1}u_{i}}{\partial
t}=\left[-\frac{\partial
C_{ni}P}{\partial\zeta_{n}}-\frac{\partial\left[C_{nj}u_{j}\right]u_{i}}{\partial\zeta_{n}}+\sqrt{\frac{Pr}{Ra}}\left(\frac{\partial}{\partial\zeta_{n}}G_{nj}\frac{\partial
u_{i}}{\partial\zeta_{j}}\right)_{n\neq j}+gT\delta_{in}\right]_{ex}\\\
+\left[\sqrt{\frac{Pr}{Ra}}\left(\frac{\partial}{\partial\zeta_{n}}G_{nj}\frac{\partial
u_{i}}{\partial\zeta_{j}}\right)_{n=j}\right]_{im}\end{split}$ (21)

$\displaystyle\begin{split}J^{-1}u_{i}^{\star}=J^{-1}u_{i}^{n}+\beta^{(1)}h^{(1)}\square+\sqrt{\frac{Pr}{Ra}}\frac{h^{(1)}}{2}\left(\frac{\partial}{\partial\zeta}\left[G_{11}\frac{\partial
u_{i}^{n}}{\partial\zeta}\right]+\frac{\partial}{\partial\zeta}\left[G_{11}\frac{\partial
u_{i}^{\star}}{\partial\zeta}\right]\right),\\\
+\sqrt{\frac{Pr}{Ra}}h^{(1)}\frac{\partial}{\partial\eta}\left[G_{22}\frac{\partial
u_{i}^{n}}{\partial\eta}\right]+\sqrt{\frac{Pr}{Ra}}h^{(1)}\frac{\partial}{\partial\xi}\left[G_{33}\frac{\partial
u_{i}^{n}}{\partial\xi}\right],\end{split}$ (22)

Here, $\square$ represents all the terms to be marched explicitly.

$\displaystyle\begin{split}J^{-1}u_{i}^{\star\star}=J^{-1}u_{i}^{\star}-\sqrt{\frac{Pr}{Ra}}\frac{h^{(1)}}{2}\frac{\partial}{\partial\eta}\left[G_{22}\frac{\partial
u_{i}^{n}}{\partial\eta}\right]+\sqrt{\frac{Pr}{Ra}}\frac{h^{(1)}}{2}\frac{\partial}{\partial\eta}\left[G_{22}\frac{\partial
u_{i}^{\star\star}}{\partial\eta}\right]\end{split}$ (23)

$\displaystyle\begin{split}J^{-1}u_{i}^{\star\star\star}=J^{-1}u_{i}^{\star\star}-\sqrt{\frac{Pr}{Ra}}\frac{h^{(1)}}{2}\frac{\partial}{\partial\xi}\left[G_{33}\frac{\partial
u_{i}^{n}}{\partial\xi}\right]+\sqrt{\frac{Pr}{Ra}}\frac{h^{(1)}}{2}\frac{\partial}{\partial\xi}\left[G_{33}\frac{\partial
u_{i}^{\star\star\star}}{\partial\xi}\right]\end{split}$ (24)

The intermediate velocity fields $u^{\star}_{i}$, $u^{\star\star}_{i}$,
$u^{\star\star\star}_{i}$ are obtained by solving a set of tridiagonal
matrices that result from the spatial discretization of the equations 22, 23,
and 24 in the $\zeta$, $\eta$, and $\xi$ directions respectively. The
intermediate velocity $u^{\star\star\star}_{i}$ is the first step in the
fractional-step scheme. We employ the Thomas Algorithm with Pipelining, as
described in the next section, to solve the tridiagonal system 22-24 to obtain
$u^{\star}$, $u^{\star\star}$, and $u^{\star\star\star}$.

###  3.5 Thomas Algorithm

Let us consider solving $A\psi=g$ for $\psi$, which is the outcome of the
spatial discretization of, for instance, 22. Here $A$ is a tridiagonal matrix
given as

$\left[\begin{array}[]{ccccccccc}b_{0~{}~{}}&c_{0~{}~{}}&&&&&&&\\\
a_{1~{}~{}}&b_{1~{}~{}}&c_{1~{}~{}}&&&&&&\\\
&a_{2~{}~{}}&b_{2~{}~{}}&c_{2~{}~{}}&&&&&\\\
&&\cdot}{}{}{}{&\cdot}{}{}{}{&\cdot}{}{}{}{&&&&\\\
&&&\cdot}{}{}{}{&\cdot}{}{}{}{&\cdot}{}{}{}{&&&\\\
&&&&a_{n-1}&b_{n-1}&c_{n-1}&&\\\ &&&&&a_{n~{}~{}}&b_{n~{}~{}}\\\
\end{array}\right]\left[{\begin{array}[]{c}\psi^{~{}}_{0~{}~{}}\\\
\psi^{~{}}_{1~{}~{}}\\\ \psi^{~{}}_{2~{}~{}}\\\ ~{}\cdot_{~{}~{}~{}}\\\
~{}\cdot_{~{}~{}~{}}\\\ \psi^{~{}}_{n-1}\\\ \psi^{~{}}_{n~{}~{}}\\\
\end{array}}\right]=\left[{\begin{array}[]{c}g^{~{}}_{0~{}~{}}\\\
g^{~{}}_{1~{}~{}}\\\ g^{~{}}_{2~{}~{}}\\\ ~{}\cdot_{~{}~{}~{}}\\\
~{}\cdot_{~{}~{}~{}}\\\ g^{~{}}_{n-1}\\\ g^{~{}}_{n~{}~{}}\\\
\end{array}}\right].$ (25)

The first two relations in 25 are,

$\displaystyle b_{0}\psi_{0}+c_{0}\psi_{1}$ $\displaystyle=g_{0}$ (26)
$\displaystyle a_{1}\psi_{0}+b_{1}\psi_{1}+c_{1}\psi_{2}$
$\displaystyle=g_{1}.$ (27)

Substituting $\psi_{0}$ from 26 into $\psi_{0}$ in 27 gives

$b^{{}^{\prime}}_{1}\psi_{1}+c_{1}\psi_{2}=g^{{}^{\prime}}_{1}$ (28)

where $b^{{}^{\prime}}_{1}=\left[b_{1}-a_{1}b^{-1}_{0}c_{0}\right]$ and
$g^{{}^{\prime}}_{1}=\left[g_{1}-a_{1}b^{-1}_{0}g_{0}\right]$. The algorithm
involves two stages, forward sweeping and backward substitution. The sub-
diagonal elements $a_{1}-a_{n}$ are removed using Gaussian elimination during
the forward sweeping step. Therefore, equation 25 takes the form:

$\left[\begin{array}[]{ccccccccc}b_{0~{}~{}}&c_{0~{}~{}}&&&&&&&\\\
0&b^{{}^{\prime}}_{1~{}~{}}&c_{1~{}~{}}&&&&&&\\\
&0&b^{{}^{\prime}}_{2~{}~{}}&c_{2~{}~{}}&&&&&\\\
&&\cdot}{}{}{}{&\cdot}{}{}{}{&\cdot}{}{}{}{&&&&\\\
&&&\cdot}{}{}{}{&\cdot}{}{}{}{&\cdot}{}{}{}{&&&\\\
&&&&0&b^{{}^{\prime}}_{n-1}&c_{n-1}&&\\\ &&&&&0&b^{{}^{\prime}}_{n~{}~{}}\\\
\end{array}\right]\left[{\begin{array}[]{c}\psi^{~{}}_{0~{}~{}}\\\
\psi^{~{}}_{1~{}~{}}\\\ \psi^{~{}}_{2~{}~{}}\\\ ~{}\cdot_{~{}~{}~{}}\\\
~{}\cdot_{~{}~{}~{}}\\\ \psi^{~{}}_{n-1}\\\ \psi^{~{}}_{n~{}~{}}\\\
\end{array}}\right]=\left[{\begin{array}[]{c}g^{~{}}_{0~{}~{}}\\\
g^{{}^{\prime}}_{1~{}~{}}\\\ g^{{}^{\prime}}_{2~{}~{}}\\\
~{}\cdot_{~{}~{}~{}}\\\ ~{}\cdot_{~{}~{}~{}}\\\ g^{{}^{\prime}}_{n-1}\\\
g^{{}^{\prime}}_{n~{}~{}}\\\ \end{array}}\right].$ (29)

At the end of the forward sweep, we can solve for
$\psi_{n}=g^{{}^{\prime}}_{n}/b^{{}^{\prime}}_{n}$ in 29. Subsequently, we
solve for $\psi_{n-1}$-$\psi_{0}$ (equations $n-1$ until $0$) as in
$\psi_{i}=(g^{{}^{\prime}}_{i}-c_{i}\psi_{i+1})/b^{{}^{\prime}}_{i}$. For the
grid in Section 3.7, periodic boundary conditions are enforced in the $\xi$
direction. Therefore, the above-mentioned Thomas algorithm is modified as
follows. Consider the discretized system $A\psi=g$ with periodicity as in 30
where $\psi_{1}=\psi_{n-1}$ and $\psi_{2}=\psi_{n}$.

$\left[\begin{array}[]{cccccccc}b_{1~{}~{}}&c_{1~{}~{}}&&&&&a_{1~{}~{}}\\\
a_{2~{}~{}}&b_{2~{}~{}}&c_{2~{}~{}}&&&&\\\
&\cdot}{}{}{}{&\cdot}{}{}{}{&\cdot}{}{}{}{&&&\\\
&&\cdot}{}{}{}{&\cdot}{}{}{}{&\cdot}{}{}{}{&&\\\
&&&\cdot}{}{}{}{&\cdot}{}{}{}{&\cdot}{}{}{}{&\\\
&&&&a_{n-1}&b_{n-1}&c_{n-1}\\\ c_{n~{}~{}}&&&&&a_{n~{}~{}}&b_{n~{}~{}}\\\
\end{array}\right]\left[{\begin{array}[]{c}\psi^{~{}}_{1~{}~{}}\\\
\psi^{~{}}_{2~{}~{}}\\\ ~{}\cdot_{~{}~{}~{}}\\\ ~{}\cdot_{~{}~{}~{}}\\\
~{}\cdot_{~{}~{}~{}}\\\ \psi^{~{}}_{n-1}\\\ \psi^{~{}}_{n~{}~{}}\\\
\end{array}}\right]=\left[{\begin{array}[]{c}g^{~{}}_{1~{}~{}}\\\
g^{~{}}_{2~{}~{}}\\\ ~{}\cdot_{~{}~{}~{}}\\\ ~{}\cdot_{~{}~{}~{}}\\\
~{}\cdot_{~{}~{}~{}}\\\ g^{~{}}_{n-1}\\\ g^{~{}}_{n~{}~{}}\\\
\end{array}}\right]$ (30)

The first step includes separating equation 30 into a tridiagonal system 31
with an additional equation 32.

$\left[\begin{array}[]{ccccccc}b_{1~{}~{}}&c_{1~{}~{}}&&&&\\\
a_{2~{}~{}}&b_{2~{}~{}}&c_{2~{}~{}}&&&\\\
&\cdot}{}{}{}{&\cdot}{}{}{}{&\cdot}{}{}{}{&&\\\
&&\cdot}{}{}{}{&\cdot}{}{}{}{&\cdot}{}{}{}{&\\\
&&&\cdot}{}{}{}{&\cdot}{}{}{}{&\cdot}{}{}{}{\\\ &&&&a_{n-1}&b_{n-1}\\\
\end{array}\right]\left[{\begin{array}[]{c}\psi^{~{}}_{1~{}~{}}\\\
\psi^{~{}}_{2~{}~{}}\\\ ~{}\cdot_{~{}~{}~{}}\\\ ~{}\cdot_{~{}~{}~{}}\\\
~{}\cdot_{~{}~{}~{}}\\\ \psi^{~{}}_{n-1}\\\
\end{array}}\right]=\left[{\begin{array}[]{c}g^{~{}}_{1~{}~{}}\\\
g^{~{}}_{2~{}~{}}\\\ ~{}\cdot_{~{}~{}~{}}\\\ ~{}\cdot_{~{}~{}~{}}\\\
~{}\cdot_{~{}~{}~{}}\\\ g^{~{}}_{n-1}\\\
\end{array}}\right]+\left[{\begin{array}[]{c}-a^{~{}}_{1~{}~{}}\\\
~{}0^{~{}}_{~{}~{}~{}}\\\ ~{}0^{~{}}_{~{}~{}~{}}\\\ ~{}\cdot_{~{}~{}~{}}\\\
~{}\cdot_{~{}~{}~{}}\\\ -c^{~{}}_{n-1}\\\ \end{array}}\right]\psi_{n}$ (31)

$c_{n}\psi_{1}+a_{n}\psi_{n-1}+b_{n}\psi_{n}=g_{n}.$ (32)

The tridiagonal matrix on the left-hand side is defined as $[A_{c}]$ and $[g]$
as the g-column matrix on the right-hand side. Let

$[\psi]=[\psi 1]+[\psi 2]\psi_{n}$ (33)

be the solution of the system 31 where

$\displaystyle[\psi 1]$ $\displaystyle=[A_{c}]^{-1}[g]$ (34)
$\displaystyle[\psi 2]$ $\displaystyle=[A_{c}]^{-1}[-a_{1}0\cdot\cdot\cdot-
c_{n-1}]^{T}.$ (35)

Substituting $\psi$ in 33 into $\psi_{1}$ and $\psi_{n-1}$ in 32 gives

$c_{n}(\psi 1_{1}+\psi 2_{1}\psi_{n})+a_{n}(\psi 1_{n-1}+\psi
2_{n-1}\psi_{n})+b_{n}\psi_{n}=g_{n}.$ (36)

Rearrange 36 for $\psi_{n}$

$\psi_{n}=\dfrac{g_{n}-c_{n}\psi 1_{1}-a_{n}\psi 1_{n-1}}{b_{n}+c_{n}\psi
2_{1}+a_{n}\psi 2_{n-1}}.$ (37)

In summary, to solve the system 30, we employ the following steps:

  1. 1.

Construct 34 and 35

  2. 2.

Solve 34 and 35 for $[\psi 1]$ and $[\psi 2]$ from index 1 to index $n-1$

  3. 3.

Substitute $[\psi 1]$ and $[\psi 2]$ into 37 and solve for $\psi_{n}$

  4. 4.

Calculate $[\psi]$ from 33 using $[\psi 1]$, $[\psi 2]$, and $\psi_{n}$

####  3.5.1 Parallel algorithm

$\begin{array}[]{c}\phantom{a}\\\ \phantom{a}\\\ \phantom{a}\\\ {\color[rgb]{0,0,1}\rotatebox{90.0}{$\text{CPU}^{i}$}}\\\ \phantom{a}\\\ \phantom{a}\\\ {\color[rgb]{1,0,0}\rotatebox{90.0}{$\text{CPU}^{i+1}$}}\\\ \phantom{a}\\\ \phantom{a}\\\ \end{array}\left|\begin{array}[]{ccccccccccc | c}{\color[rgb]{0.5,0.5,0}b^{{}^{\prime}i-1}_{s~{}~{}}}&{\color[rgb]{0.5,0.5,0}c^{i-1}_{s~{}~{}}}&&&&&&&&&&\cdot}{}{}{}{\\\ \hline\cr 0&b^{{}^{\prime}i~{}~{}}_{s~{}~{}}&c^{i~{}~{}}_{s~{}~{}}&&&&&&&&&g^{{}^{\prime}i~{}~{}}_{s~{}~{}}\\\ &0&b^{{}^{\prime}i~{}~{}}_{s+1}&c^{i~{}~{}}_{s+1}&&&&&&&&g^{{}^{\prime}i~{}~{}}_{s+1}\\\ &&\cdot}{}{}{}{&\cdot}{}{}{}{&\cdot}{}{}{}{&&&&&&&\cdot}{}{}{}{\\\ &&&0&{\color[rgb]{0,0,1}b^{{}^{\prime}i~{}~{}}_{e~{}~{}}}&{\color[rgb]{0,0,1}c^{i~{}~{}}_{e~{}~{}}}&&&&&&{\color[rgb]{0,0,1}g^{{}^{\prime}i~{}~{}}_{e~{}~{}}}\\\ \hline\cr&&&&a^{i+1}_{s~{}~{}}&b^{i+1}_{s~{}~{}}&c^{i+1}_{s~{}~{}}&&&&&{\color[rgb]{1,0,0}g^{i+1}_{s~{}~{}}}\\\ &&&&&a^{i+1}_{s+1}&b^{i+1}_{s+1}&c^{i+1}_{s+1}&&&&g^{i+1}_{s+1}\\\ &&&&&&\cdot}{}{}{}{&\cdot}{}{}{}{&\cdot}{}{}{}{&&&\cdot}{}{}{}{\\\ &&&&&&&a^{i+1}_{e~{}~{}}&{\color[rgb]{1,0,0}b^{i+1}_{e~{}~{}}}&{\color[rgb]{1,0,0}c^{i+1}_{e~{}~{}}}&&{\color[rgb]{1,0,0}g^{i+1}_{e~{}~{}}}\\\ \hline\cr&&&&&&&&{\color[rgb]{1,0,0}a^{i+1}_{s~{}~{}}}&\cdot}{}{}{}{&\cdot}{}{}{}{&\cdot}{}{}{}{\\\ \end{array}\right|$ (38) $\begin{array}[]{c}\phantom{a}\\\ \phantom{a}\\\ \phantom{a}\\\ {\color[rgb]{0,0,1}\rotatebox{90.0}{$\text{CPU}^{i}$}}\\\ \phantom{a}\\\ \phantom{a}\\\ {\color[rgb]{1,0,0}\rotatebox{90.0}{$\text{CPU}^{i+1}$}}\\\ \phantom{a}\\\ \phantom{a}\\\ \end{array}\left|\begin{array}[]{ccccccccccc | c}b^{{}^{\prime}i-1}_{s~{}~{}}&c^{i-1}_{s~{}~{}}&&&&&&&&&&\cdot}{}{}{}{\\\ \hline\cr&b^{{}^{\prime}i~{}~{}}_{s~{}~{}}&c^{i~{}~{}}_{s~{}~{}}&&&&&&&&&g^{{}^{\prime}i~{}~{}}_{s~{}~{}}\\\ &&b^{{}^{\prime}i~{}~{}}_{s+1}&c^{i~{}~{}}_{s+1}&&&&&&&&g^{{}^{\prime}i~{}~{}}_{s+1}\\\ &&&\cdot}{}{}{}{&\cdot}{}{}{}{&&&&&&&\cdot}{}{}{}{\\\ &&&&b^{{}^{\prime}i~{}~{}}_{e~{}~{}}&c^{i~{}~{}}_{e~{}~{}}&&&&&&g^{{}^{\prime}i~{}~{}}_{e~{}~{}}\\\ \hline\cr&&&&&b^{{}^{\prime}i+1}_{s~{}~{}}&c^{i+1}_{s~{}~{}}&&&&&{\color[rgb]{1,0,0}\psi^{i+1}_{s~{}~{}}}\\\ &&&&&&b^{{}^{\prime}i+1}_{s+1}&c^{i+1}_{s+1}&&&&\psi^{i+1}_{s+1}\\\ &&&&&&&\cdot}{}{}{}{&\cdot}{}{}{}{&&&\cdot}{}{}{}{\\\ &&&&&&&&b^{{}^{\prime}i+1}_{e~{}~{}}&c^{i+1}_{e~{}~{}}&&\psi^{i+1}_{e~{}~{}}\\\ \hline\cr&&&&&&&&&\cdot}{}{}{}{&\cdot}{}{}{}{\\\ \end{array}\right|$ (39)

Fig. 6:  Example of a 4$\times$4$\times$4 CPU topology in a single
computational domain. Numbers in the figure represent an individual CPU’s rank
(id).

Figure 6 demonstrates an example 4$\times$4$\times$4 central processing unit
(CPU) topology in a single computational domain. The spatial discretization of
22, 23, or 24 over several CPUs results in tridiagonal matrices, and vectors.
The tridiagonal matrices constructed from the discretization of 24 span over
the entire $\xi$ space index, e.g. over CPU0-3 illustrated by the solid red
line in the figure . Similarly, the tridiagonal matrices constructed from the
discretization of 23 span over the entire $\eta$ space index (e.g. over
CPU19-31), as indicated by a solid blue line. Considering 38, forward sweeping
starts at CPU0. Once the sweeping reaches the interface between CPUi and
CPUi+1, CPUi sends $b^{{}^{\prime}i}_{e}$, $c^{i}_{e}$, and
$g^{{}^{\prime}i}_{e}$ to CPUi+1. Then, CPUi+1 continues to carry out the
sweeping by sending data $b^{{}^{\prime}i+1}_{e}$, $c^{i+1}_{e}$, and
$g^{{}^{\prime}i+1}_{e}$ to CPUi+2 and so on. After the forward sweeping is
finalized, backward substitution starts, and reverse sweeping is performed, as
shown in 39. However, the only information being sent from CPUi+1 to CPUi is
$\psi^{i+1}_{s}$.

For a periodic system, we use the following steps:

  1. 1.

Construct 34 and 35

  2. 2.

Use the parallel Thomas subroutine to solve 34 and 35 for $[\psi 1]$ and
$[\psi 2]$

  3. 3.

CPU0 owning the first block (contains node 1), sends $\psi 1_{1}$ and $\psi
2_{1}$ to CPUN that owns the last block (contains node n)

  4. 4.

CPUN calculates $\psi_{n}$ and broadcasts $\psi_{n}$ to every CPU that owns a
subsystem of 30

  5. 5.

Every CPU calculates $[\psi]$ from 33

Notice that by splitting 30 into 31 and 32, CPUN solves the tridiagonal system
33 which has size one element less than the others.

####  3.5.2 Pipelining

In the previous subsection, we summarize how to solve a tridiagonal system in
parallel. It is done simply by completing the forward/backward sweep and
sending data to the proper neighbor in order to continue marching. CPUi that
finishes the forward sweep sends data to CPUi+1 until the last block is
reached. Generally, each CPU can be responsible for thousands of tridiagonal
subsystems contained in a single subdomain (or ‘a block’). This subsection
summarizes how to solve such a big system efficiently.

Consider a computational domain containing $(n_{\zeta},n_{\eta},n_{\xi})$ grid
points. For the sake of simplicity, the domain is equally decomposed only in
the $\eta$ direction into $NJ$ blocks so that CPU0 occupies block $0$, CPU1
occupies block $1$, and so on. Thus, each CPU owns a block of size
$(n_{\zeta},n_{\eta}/N_{\eta},n_{\xi})$; given that $n_{\eta}/N_{\eta}$ is, by
design, an integer. Supposing that we choose to perform an implicit marching
in the $\eta$ direction, the resulting tridiagonal matrix is subdivided into
$N_{\eta}$ sections. The easiest, though the least efficient, way to solve
these systems is to let CPU0 solve all of its subsystems across the
$(n_{\zeta},n_{\eta}/N_{\eta},n_{\xi})$ grid before sending data to CPU1. That
is, CPU0 performs a forward sweep at cell $(1,1\rightarrow
n_{\eta}/N_{\eta},1)$, at cell $(1,1\rightarrow n_{\eta}/N_{\eta},2)$, and so
on until cell $(n_{\zeta},1\rightarrow n_{\eta}/N_{\eta},n_{\xi})$. Next, CPU0
packs the plane data with $n_{\zeta}*n_{\xi}*3$ elements (recall
$b^{{}^{\prime}i}_{e}$, $c^{i}_{e}$, and $g^{{}^{\prime}i}_{e}$ in the
previous subsection) at $(1:n_{\zeta},n_{\eta}/N_{\eta},1:n_{\xi})$ and sends
it to CPU1. Following the same process for the subsequent CPUs until
CPU${}^{N_{\eta}-1}$ is reached, the backward substitution is carried out in
the same way from CPU${}^{N_{\eta}-1}$ to CPU0. The obvious drawback is that
only one CPU operates at a given time, and the whole process will be even
slower than the serial version since there is additional communication
overhead.

Fig. 7: Illustration of pipelining with the parallel Thomas algorithm
following Figure VI.16 of Taylor [27]. Forward sweep is achieved in steps
(a-d), while steps (e-g) depict back-substitution.

Pipelining is employed in an attempt to minimize the number of idle CPUs while
optimizing communication overhead. In essence, rather than sweeping across the
$(n_{\zeta},n_{\eta}/N_{\eta},n_{\xi})$ grid all at once, each CPU performs
the sweeps only for a portion of the grid and shares data with its
neighbouring CPU in the sweep direction downstream in a forward sweep and
upstream for a backward sweep. A portion of the grid can be chosen for the
first CPU, with the others obeying the same portion. We give an example of
pencil-type pipelining. Consider figure 7 and the following steps:

  1. 1.

CPU0, process rank 0 in the figure, performs the forward sweep in the $\eta$
direction at cell $(i,k)=(1,1)$ from $(i,j,k)=(1,1,1)$ to
$(i,j,k)=(1,n_{\eta}/N_{\eta},1)$; here $i$ and $k$ are dummy indices pointing
to a grid location in $\zeta$ and $\xi$ directions, respectively. CPU0 then
repeats the forward sweep until $(i,k)=(1,n_{\xi})$. Notice that the forward
sweep is in the $\eta$ direction, but the ‘pencil’ aligns in the $\xi$
direction. At this point, CPU0 packs and passes data to CPU1. The data is of
size $n_{\xi}*3$ elements containing $b^{{}^{\prime}0}_{n_{\eta}/N_{\eta}}$,
$c^{0}_{n_{\eta}/N_{\eta}}$, and $g^{{}^{\prime}0}_{n_{\eta}/N_{\eta}}$ for
each $k\in[1,n_{\xi}]$ (with 1-element width in the $\zeta$ direction, hence
the word ‘pencil’).

  2. 2.

CPU1 continues the forward sweep while CPU0 starts solving the new tridiagonal
system by shifting 1 step from the first block in the $\zeta$, which is the
‘slide’ direction. The ‘slide’ and ‘pencil’ directions can be swapped.

  3. 3.

CPU1 passes data to CPU2 for the sliding index $i=1$, receives data from CPU0
at the sliding index $i=2$, and continues the forward sweep.

  4. 4.

The same process is carried out until CPU${}^{N_{\eta}-1}$ reaches the slide
index $i=n_{\zeta}$.

  5. 5.

CPU${}^{N_{\eta}-1}$ starts the backward sweep at the sliding index
$i=n_{\zeta}$, shares data of size $n_{\xi}*1$-element containing
$\zeta^{N_{\eta}-1}_{1}$ for each $k\in[1,n_{\xi}]$ with CPU${}^{N_{\eta}-2}$,
and starts the backward sweep at the sliding index $i=n_{\zeta}-1$.

  6. 6.

The backward sweeping process is carried out in the same way as the forward
sweep.

  7. 7.

Solving the system of tridiagonal matrices is finalized after CPU0 finishes
the backward sweep at the sliding index $i=1$.

####  3.5.3 Handling shell cut

Using the parallel Thomas algorithm with pipelining, we have been able to
solve 22, 23, and 24 in parallel for $u^{\star}$, $u^{\star\star}$, and
$u^{\star\star\star}$. The grid used for the solver before it is rotated about
the $x_{1}$-axis is shown in figure 2a. The directions parallel and
perpendicular to the body surface are denoted by $\zeta$ and $\eta$,
respectively. Figure 2b represents the transformed coordinate, that is
obtained using Jacobi transformation. The top and bottom edges of the domain,
parallel to $x_{1}$-axis in figure 2a, are indicated by the phrase ’Branch
cut’ (AB and CD). They correspond to the left and right sides of the
transformed domain, which are parallel to the $\eta$ direction, indicated by
the word ’shell cut’. The body surface seen in the curvilinear domain in
figure 2a is transformed to the top and bottom body surface of the transformed
domain 2b. The words ’shell cut’ or ’Branch cut’ represent a shared interface
among CPUs that cuts through the centerline. A tridiagonal system in the
$\zeta$ direction created by discretizing 22 is interrupted by the shell cut
on both the left and right sides. To handle the shell cut, the tridiagonal
system from one side of the cut is merged with the system on the opposite
side. Therefore, the resulting system is twice as large as the system without
the cut.

Fig. 8:  Solving a tridiagonal system across a shell cut. Here (A, A’) and (B,
B’) are the two pairs of surfaces across the left and right shell cuts,
respectively, as per the arrangement of processors shown in figure 6.

For example, consider solving a system in the $\zeta$ direction in figure 8.
In the forward sweeping, CPU1 with point A starts solving from this point A
and then passes the information to CPU17 until the forward sweeping reaches
CPU49. Then, CPU49, with point B passes the information to CPU51 that has
point B${}^{{}^{\prime}}$ on the other side of the cut. CPU51 keeps performing
the forward sweep by sending data to CPU35 and so on until the sweep reaches
$A^{{}^{\prime}}$ owned by CPU3. The backward substitution follows the same
procedure by starting from point A${}^{{}^{\prime}}$ and marching until the
substitution reaches back to point A. Forward and backward sweeping are done
using the pencil-type pipeline Thomas algorithm explained previously.

###  3.6 Pressure Correction

The generalized curvilinear solver uses a combination of the ADI-CN-RKW3
methods to obtain $u^{\star\star\star}$, the intermediate velocity. The
remaining procedure in the fractional step method is to remove the divergence
residual from the projected velocity $u^{\star\star\star}$ at the end of each
sub-RKW3 step (denoted as PC in figure 5). This step requires correcting the
pressure to account for the divergenceless field. Rewriting equation 6 as 40;
where $\Box$ represents the advection, the diffusion, and the baroclinic
terms. 40 is temporally discretized into 41 and 42. Here,
$u^{\star\star\star}_{i}$ denotes velocity at the third step of the ADI, and
${h}$ is a sub-time step of RKW3.

$\frac{\partial J^{-1}u_{i}}{\partial t}=\Box-\frac{\partial
C_{ji}P}{\partial\zeta_{j}}$ (40)
$\frac{J^{-1}u^{\star\star\star}_{i}-J^{-1}u^{n}_{i}}{{h}}=\Box^{n}-\frac{\partial
C_{ji}P^{n}}{\partial\zeta_{j}}$ (41)
$\frac{J^{-1}u^{n+{h}}_{i}-J^{-1}u^{n}_{i}}{{h}}=\Box^{n}-\frac{\partial
C_{ji}P^{n+{h}}}{\partial\zeta_{j}}$ (42)

42-41 gives:

$\frac{J^{-1}u^{n+{h}}_{i}-J^{-1}u^{\star\star\star}_{i}}{{h}}=-\frac{\partial
C_{ji}\delta P^{{h}}}{\partial\zeta_{j}}$ (43)

Taking divergence of 43 gives 44. Note that $\partial_{i}u^{n+{h}}_{i}=0$.

$\frac{1}{{h}}\frac{\partial}{\partial\zeta_{j}}\left[\frac{\partial\zeta_{j}}{\partial
x_{i}}J^{-1}u^{\star\star\star}_{i}\right]=\frac{\partial}{\partial\zeta_{k}}\frac{\partial\zeta_{k}}{\partial
x_{i}}\left[\frac{\partial}{\partial\zeta_{j}}J^{-1}\frac{\partial\zeta_{j}}{\partial
x_{i}}\delta P^{{h}}\right]$ (44)

This yields the Poisson equation 45 for pressure correction $\delta P^{{h}}$

$\frac{\partial}{\partial\zeta_{i}}\frac{\partial}{\partial\zeta_{j}}\left[G_{ij}{h}\delta
P^{{h}}\right]=\frac{\partial}{\partial\zeta_{j}}\left[C_{ji}u^{\star\star\star}_{i}\right]$
(45)

Removing divergence from the $u^{\star\star\star}$ field is done by solving 45
for $\delta P^{\bar{h}}$ and computing $u^{n+\bar{h}}_{i}$ using 46.

$u^{n+\bar{h}}_{i}=u^{\star\star\star}_{i}-\frac{1}{J^{-1}}\frac{\partial}{\partial\zeta_{j}}\left[C_{ji}\bar{h}P^{\bar{h}}\right]$
(46)

The Poisson equation for pressure correction is solved using the Semi-
Coarsening Multigrid routine in the HYPRE library [28]. The divergence-free
field $u^{n+\bar{h^{(1)}}}$ marks the end of RKW3 first sub step. We follow
the same procedure until $u^{n+\bar{h^{(1)}}+\bar{h^{(2)}}+\bar{h^{(3)}}}$ =
$u^{n+1}$ is obtained. Figure 5 illustrates the entire process. HYPRE is a
library of scalable linear solvers and multigrid methods [29] and [30].
Generalized curvilinear solver utilizes two solvers provided by HYPRE: 1) SMG,
a parallel semi-coarsening multigrid solver for linear systems [28] and 2)
BoomerAMG, a parallel implementation of the algebraic multigrid method [31].

###  3.7 Simulation details

We have performed simulations for three Rayleigh numbers, $Ra=10^{5},10^{7}$
and $10^{8}$, keeping the Prandtl number constant at $Pr=1$. The radius ratio
is also kept to be constant at $\Gamma=0.6$, and an inverse square-law profile
of gravity with the radius $g=(r_{o}/r)^{2}$ is assumed. These choices of
parameters are aimed to facilitate comparison with Gastine et al. [14]. The
number of gridpoints used in each direction for different $Ra$ is given in
table 2. The physical curvilinear grid is clustered in the radial direction
near the boundaries to resolve the boundary layers near the solid surfaces
before performing the Jacobi transformation. The clustering function is given
below.

$r(j)=\frac{\tanh\left[rx_{2}\left(\frac{j-1}{nx_{2}}-\frac{1}{2}\right)\right]}{2\tanh\left(\frac{rx_{2}}{2}\right)},$
(47)

here, $rx_{2}$ is the stretching factor, and $nx_{2}$ is the number of grid
divisions in the radial direction.

After the transformation, all the grid spacings are unity, and the information
about grid stretching is provided effectively through the elongation matrix.
At the bottom and top surfaces, a no-slip boundary condition is used
($u_{1}=u_{2}=u_{3}=0$), while the temperatures are fixed at the bottom
($T=1$) and top ($T=0$) surfaces to impose an unstable gradient for
maintaining thermal convection. The periodic boundary condition is used for
all the variables in the $\xi$ directions. The ”shell-cut” boundary conditions
are used in the $\zeta$ direction, as explained before in section 3.5.3. All
simulations are started with $u_{i}=0$ and small random perturbations in the
temperature field.

We use the following notations for the surface, volume, and time-averaged
quantities.

$\left\langle
f\right\rangle_{s}=\frac{1}{4\pi}\int_{0}^{2\pi}\int_{0}^{\pi}f\;\sin{\theta}\;d\theta\;d\phi,$
(48)

$\left\langle
f\right\rangle=\frac{1}{V}\int_{r_{i}}^{r_{o}}\int_{0}^{2\pi}\int_{0}^{\pi}f\;r^{2}\;\sin{\theta}\;d\theta\;d\phi\;dr,$
(49)

$\overline{f}=\frac{1}{\tau}\int_{t_{0}}^{t_{0}+\tau}f\;d\tau,$ (50)

where $V=\frac{4}{3}\pi(r_{o}^{3}-r_{i}^{3})$. All the statistical quantities
are averaged in time for at least $100$ free fall time ($t_{f}=d/u_{f}$) units
after the simulation reaches a steady state. Heat transport in the spherical
shell is quantified by the Nusselt number $Nu$, which is defined as

$Nu=\frac{\overline{\left\langle
u_{r}T\right\rangle_{s}}-\frac{1}{\sqrt{RaPr}}\frac{\mathrm{d}\rho}{\mathrm{d}r}}{-\frac{1}{\sqrt{RaPr}}\frac{\mathrm{d}T_{c}}{\mathrm{dr}}}=-\Gamma\frac{\mathrm{d}\rho}{\mathrm{d}r}\left(r=r_{i}\right)=-\frac{1}{\Gamma}\frac{\mathrm{d}\rho}{\mathrm{d}r}\left(r=r_{o}\right),$
(51)

where $\rho(r)=\overline{\langle T\rangle_{s}}$, $\langle\rangle_{s}$
represents the average over the spherical surface 48 and the overbar
represents the time average 50. Here, $T_{c}$ is the conductive temperature
profile for spherical shells with isothermal boundaries given by 53. The
thermal conduction equation for a spherical shell with isothermal boundary
condition is given by

$\frac{\mathrm{d}}{\mathrm{d}r}\left(r^{2}\frac{\mathrm{d}T_{c}}{\mathrm{~{}d}r}\right)=0,\quad
T_{c}\left(r=r_{i}\right)=1,\quad T_{c}\left(r=r_{o}\right)=0,$ (52)

which yields

$T_{c}(r)=\frac{\Gamma}{(1-\Gamma)^{2}}\frac{1}{r}-\frac{\Gamma}{1-\Gamma}.$
(53)

(a)

(b)

Fig. 9: Comparison of surface and time-averaged non-dimensional radial
temperature profile for the cases at (a) $Ra=10^{7}$ and (b) $Ra=10^{8}$.

##  4 Results

This section summarizes the results of the simulations listed in table 2. We
validate our results with those of [14], and further demonstrate the closure
of the turbulent kinetic energy budget.

###  4.1 Validation

$Ra$ | Grid | $Nu$ | $Nu$ | $\delta^{T}_{i}/\delta^{T}_{o}$ | $\delta^{u}_{i}/\delta^{u}_{o}$
---|---|---|---|---|---
| ($n_{\zeta}\times n_{\eta}\times n_{\xi}$) | (present DNS) | (Gastine et al. [14]) |  |
$10^{5}$ | $256\times 128\times 256$ | 4.60 | 4.71 | 0.095/0.132 | 0.139/0.209
$10^{7}$ | $512\times 256\times 512$ | 17.47 | 17.07 | 0.018/0.031 | 0.076/0.102
$10^{8}$ | $512\times 256\times 512$ | 36.50 | 33.54 | 0.009/0.015 | 0.063/0.084
Table 2: Summary table for different $Ra$ along with the grid used for each
case. Here, Prandtl number $Pr=1$, gravity profile, $g=(r_{o}/r)^{2}$, and
radius ratio $\Gamma=0.6$ has been used for all the runs.

We compare the non-dimensionalized temperature profile variation in the radial
direction in figure 9. The radial temperature variation is found to be
matching with that of Gastine et al. [14]. Additionally, we compare the $Nu$
obtained from our solver with the values reported by Gastine et al. [14] at
the same $Ra$. For the calculation of $Nu$, we use equation 51 with $r_{i}=3$
and $r_{o}=5$. From the table 2, it can be observed that $Nu$ matches very
well with the values from the Gastine et al. [14].

###  4.2 Flow visualization

A comparison of the qualitative features of the instantaneous flow and the
thermal field with an increase of $Ra$ from $10^{5}$ to $10^{8}$ is presented
in figure 10. For lower thermal forcing at $Ra=10^{5}$, the plumes generated
from the boundaries span the radial extent of the domain, as seen from 10(a).
However, $Ra=10^{8}$, as shown in 10(b), the plumes are much smaller with a
well-mixed interior. With the increase in $Ra$, turbulence increases,
accompanied by higher mixing and generation of smaller scales. The higher $Ra$
case will have a negligible temperature gradient in bulk due to enhanced
mixing. In figure 10(c), the alternating regions with positive and negative
radial velocities indicate the presence of structures similar to convective
rolls, while figure 10(d) exhibit the presence of small-scale plumes near the
boundary at higher $Ra$.

(a)

(b)

(c)

(d)

Fig. 10: Instantaneous snapshots of (a,b) the thermal field and (c,d) radial
velocity field observed at Rayleigh numbers (a,c) $Ra=10^{5}$ and (b,d)
$Ra=10^{8}$.

###  4.3 Boundary layer asymmetry

The thermal boundary layer thicknesses ($\delta^{T}_{i}$, inner
$\delta^{T}_{o}$, outer) are defined as the distance of the local maximums in
the $T_{rms}$ 54 profile from the inner and the outer walls, respectively. The
velocity boundary layer thicknesses ($\delta^{u}_{i}$, inner $\delta^{u}_{o}$,
outer) are evaluated similarly from the location of the local maximums in the
horizontal velocity profile, $u_{h}$, 55 [32]. From the values of boundary
layer thickness at the inner and outer boundary from the table 2, it is
visible that there is an asymmetry in the temperature profile at both the
inner and outer radius. The total heat flowing in through the inner surface
should flow out from the outer surface for thermal equilibrium. In conjunction
with the inner spherical shell area being less than the outer spherical shell
area, the temperature drop is higher at the inner boundary than at the outer
boundary 14. It is also visible that with an increase in $Ra$, steepening of
the temperature profile near the boundaries occurs, as shown in the insets of
11.

$T_{rms}(r)=\sqrt{\left\langle(T-\langle T\rangle_{s})^{2}\right\rangle_{s}}$
(54) $u_{h}(r)=\sqrt{\left\langle
u_{\theta}^{2}+u_{\phi}^{2}\;\right\rangle_{s}}$ (55)

Fig. 11: Temperature variations observed across different $Ra$. An enlarged
view in the inset to demonstrate the increase in gradient with the increase in
$Ra$).

###  4.4 Turbulent kinetic energy budget

We discuss the turbulent kinetic energy ($t.k.e.$) budget in RBC in this
section to review not only the kinetic energy balance but also the adequacy of
the resolution of the present simulations. The budget can be expressed as
follows:

$\left\langle\frac{\mathrm{d\mathcal{K}}}{\mathrm{d}t}\right\rangle=\left\langle\mathcal{B}\right\rangle-\left\langle\epsilon_{\nu}\right\rangle,$
(56)

where,

$\left\langle\mathcal{K}\right\rangle=\left\langle\frac{1}{2}u_{i}u_{i}\right\rangle,\;\left\langle\mathcal{B}\right\rangle=\left\langle
gu_{r}T\right\rangle,\;\left\langle\epsilon_{\nu}\right\rangle=\sqrt{\frac{Pr}{Ra}}\left\langle(\nabla\times
u)^{2}\right\rangle.$ (57)

In the RHS of equation 56, the buoyancy flux
$\left\langle\mathcal{B}\right\rangle$ is the source term that converts the
available potential energy to turbulent kinetic energy to drive the convective
motions. This $t.k.e.$ is converted to internal energy by the viscous
dissipation term $\left\langle\epsilon_{\nu}\right\rangle$,[33] which acts as
a sink. Buoyancy flux averaged over the whole spherical volume can be
expressed as,

$\left\langle\mathcal{B}\right\rangle=\frac{4\pi}{V}\int_{r_{i}}^{r_{o}}gr^{2}\overline{\left\langle
u_{r}T\right\rangle_{s}}\;dr.$ (58)

After substituting $Nu$ from 51 and $T_{c}$ from 53, we obtain,

$\left\langle\mathcal{B}\right\rangle=\frac{3}{1+\Gamma+\Gamma^{2}}\frac{1}{\sqrt{RaPr}}(Nu-1)=\sqrt{\frac{Pr}{Ra}}\;\overline{\left\langle(\nabla\times
u)^{2}\right\rangle},$ (59)
$\chi_{\epsilon_{\nu}}=\frac{\sqrt{\frac{Pr}{Ra}}\;\overline{\left\langle(\nabla\times
u)^{2}\right\rangle}}{\frac{3}{1+\Gamma+\Gamma^{2}}\frac{1}{\sqrt{RaPr}}(Nu-1)}.$
(60)

The evolution of the $t.k.e.$ budget for the case $Ra=10^{7}$ is shown in
figure 12(a). We evaluate the volume-averaged $t.k.e.$ budget terms when the
simulation becomes statistically stationary ($t\geqslant t_{0}$, where $t_{0}$
represents the starting time of the simulations). The balance term signifies
the difference between the left and right-hand sides of equation 56. This
quantity remains smaller than $5\%$ of $\langle\mathcal{B}\rangle$, indicating
sufficient resolution achieved in the simulation to dissipate all the kinetic
energy. To further quantify the spatial resolution of the numerical model, the
viscous dissipation ratio ($\chi_{\epsilon_{\nu}}$) defined by 60 is also
tested for its closeness to unity. The table 3 shows that
$\chi_{\epsilon_{\nu}}$ value is near unity for all cases, signifying good
spatial resolution. To further test the adequacy of the resolution, the radial
grid spacing is compared against the Kolmogorov scale ($l_{\eta}$) defined by,
$l_{\eta}=\left(Pr/Ra\right)^{\frac{3}{8}}\left(1/\overline{\epsilon_{\nu}}\right)^{\frac{1}{4}}$.
As seen from the figure 12(b), the radial grid spacing normalized by
$l_{\eta}$ is near unity for all the $Ra$ cases near the walls, indicating
appropriate wall resolution. As seen in figure 12(b), the normalized spacing
stays below $4$ for all the cases, which is sufficient for accurate
calculation of second-order correlations [34].

(a)

(b)

Fig. 12: (a) The $t.k.e.$ budget terms for $Ra=10^{7}$, and (b) the radial variation of grid spacing ($\Delta r$) normalized by the Kolmogorov scale ($l_{\eta}$) as estimated from the spatially averaged dissipation ($\left\langle\epsilon_{\nu}\right\rangle_{s}$). $Ra$ | $\overline{\left\langle\mathcal{B}\right\rangle}$ | $\overline{\left\langle\epsilon_{\nu}\right\rangle}$ | $\chi_{\epsilon_{\nu}}$
---|---|---|---
$10^{5}$ | $1.90\times 10^{-2}$ | $-1.69\times 10^{-2}$ | 0.97
$10^{7}$ | $9.25\times 10^{-3}$ | $-8.17\times 10^{-3}$ | 1.03
$10^{8}$ | $6.16\times 10^{-3}$ | $-5.54\times 10^{-3}$ | 1.02
Table 3: The turbulent kinetic energy budget terms and the viscous dissipation
ratio, $\chi_{\epsilon_{\nu}}$, at different $Ra$ with $Pr=1$.

##  5 Conclusion

This present investigation discusses the development of a generalized
curvilinear solver for spherical Rayleigh-Bénard convection. Using the Jacobi
transformation, the solver transforms a curvilinear domain into a Cartesian
domain, and a set of modified governing equations are solved in the Cartesian
domain. The solver uses a second-order central differencing scheme for spatial
discretization, while for temporal discretization, a combined marching scheme
of ADI-CN-RKW3 is used. A parallel Thomas algorithm with pipelining is used
for solving the tridiagonal system, which is more efficient and faster as it
reduces the idle time for CPUs. In order to remove the divergence residual
from the projected velocity in the intermediate field of the fractional step
method, Semi-coarsening Multigrid(SMG) routine from the HYPRE library for the
pressure correction is used. The solver simulates three Rayleigh number ($Ra$)
cases, namely, $10^{5}$,$10^{7}$, and $10^{8}$. The primary emphasis is given
to the ability of the solver to predict the heat transfer, quantified by the
Nusselt number, $Nu$. Comparing $Nu$ obtained from the numerical simulation
with the expected value is not a reliable criterion to assess its validity
because even the under-resolved schemes show good closeness with the $Nu$
while producing temperature fields with strong nonphysical oscillations [35].
Due to this fact, the solver is not only validated for its $Nu$ but also with
the radial temperature profiles from Gastine et al. [14]. The radial
temperature profile and $Nu$ obtained from our solver demonstrate a good match
with the results from Gastine et al. [14]. To further test the spatial
resolution, we check on the viscous dissipation ratio’s closeness to unity and
turbulent kinetic energy budget closure. The $t.k.e.$ budget reveals good
closure for all the $Ra$ cases considered. With increased $Ra$, the
temperature profile near the boundaries becomes steeper. This is because the
fluid near the boundaries is subjected to strong thermal gradients, generating
a large buoyancy force that drives the flow. Therefore, the steepening of the
temperature profile near the boundaries is evidence of the buoyancy-induced
strong convective flow with increased $Ra$. For a particular $Ra$, the
temperature profile shows asymmetry due to the difference in area between the
spherical inner and outer shell.

The majority of the computational methods developed for spherical shell
Rayleigh-Bénard convection employ spherical harmonic decomposition of the
solution variables in the angular coordinates $(\theta,\phi)$ while using
finite difference or Chebyshev polynomials in the radial direction [36, 37,
38, 39, 40, 41, 42]. A notable exception to this rule is reported in the work
by Kageyama et al. [43]. However, all these methods were developed for
perfectly spherical geometries. The novelty of the present solver lies in its
capability to account for the effects of non-spherical geometries in planetary
core convection.

Our ongoing work is focused primarily on extending the present solver to
include the effects of rotation and magnetic field. Future extensions, with
further model improvements, should reveal the possible effects of a non-
spherical geometry, not only on the convective patterns but also on the self-
generated magnetic field in global numerical dynamo simulations.

## References

  * Hartmann et al. [2001] D. L. Hartmann, L. A. Moy, Q. Fu,  Tropical convection and the energy balance at the top of the atmosphere,  J. Clim. 14 (2001) 4495–4511.
  * Roberts and King [2013] P. Roberts, E. King,  On the genesis of the earth’s magnetism,  Rep. Prog. Phys. 76 (2013) 096801.
  * Rüdiger and Hollerbach [2006] G. Rüdiger, R. Hollerbach, The magnetic universe: geophysical and astrophysical dynamo theory, John Wiley & Sons, 2006.
  * Chillà and Schumacher [2012] F. Chillà, J. Schumacher,  New perspectives in turbulent rayleigh-bénard convection,  Eur. Phys. J. E 35 (2012) 1–25.
  * Incropera [1988] F. P. Incropera,  Convection heat transfer in electronic equipment cooling (1988).
  * Incropera et al. [1996] F. P. Incropera, D. P. DeWitt, T. L. Bergman, A. S. Lavine, et al., Fundamentals of heat and mass transfer, volume 6, Wiley New York, 1996\.
  * Naskar and Pal [2022a] S. Naskar, A. Pal,  Direct numerical simulations of optimal thermal convection in rotating plane layer dynamos,  Journal of Fluid Mechanics 942 (2022a).
  * Naskar and Pal [2022b] S. Naskar, A. Pal,  Effects of kinematic and magnetic boundary conditions on the dynamics of convection-driven plane layer dynamos,  Journal of Fluid Mechanics 951 (2022b) A7.
  * Ahlers et al. [2009] G. Ahlers, S. Grossmann, D. Lohse,  Heat transfer and large scale dynamics in turbulent rayleigh-bénard convection,  Rev. Mod. Phys. 81 (2009) 503.
  * Jones [2011] C. A. Jones,  Planetary magnetic fields and fluid dynamos,  Ann. Rev. Fluid Mech. 43 (2011) 583–614.
  * Wolstencroft et al. [2009] M. Wolstencroft, J. H. Davies, D. R. Davies,  Nusselt–rayleigh number scaling for spherical shell earth mantle simulation up to a rayleigh number of $10^{9}$,  Phys. Earth Planet. Inter. 176 (2009) 132–141.
  * Gastine et al. [2016] T. Gastine, J. Wicht, J. Aubert,  Scaling regimes in spherical shell rotating convection,  J. Fluid Mech. 808 (2016) 690–732.
  * Aurnou et al. [2015] J. M. Aurnou, M. A. Calkins, J. S. Cheng, K. Julien, E. M. King, D. Nieves, K. M. Soderlund, S. Stellmach,  Rotating convective turbulence in earth and planetary cores,  Phys. Earth Planet. Inter. 246 (2015) 52–71.
  * Gastine et al. [2015] T. Gastine, J. Wicht, J. M. Aurnou,  Turbulent rayleigh-bénard convection in spherical shells,  Journal of Fluid Mechanics 778 (2015) 721–764.
  * Mound and Davies [2017] J. E. Mound, C. J. Davies,  Heat transfer in rapidly rotating convection with heterogeneous thermal boundary conditions,  J. Fluid Mech. 828 (2017) 601–629.
  * Long et al. [2020] R. S. Long, J. E. Mound, C. J. Davies, S. M. Tobias,  Scaling behaviour in spherical shell rotating convection with fixed-flux thermal boundary conditions,  Journal of Fluid Mechanics 889 (2020) A7.
  * Yadav and Bloxham [2020] R. K. Yadav, J. Bloxham,  Deep rotating convection generates the polar hexagon on saturn,  Proc. Natl. Acad. Sci. USA 117 (2020) 13991–13996.
  * Yadav et al. [2020] R. K. Yadav, M. Heimpel, J. Bloxham,  Deep convection–driven vortex formation on jupiter and saturn,  Sci. Adv. 6 (2020) eabb9298.
  * Korre and Featherstone [2021] L. Korre, N. A. Featherstone,  On the dynamics of overshooting convection in spherical shells: Effect of density stratification and rotation,  The Astrophysical Journal 923 (2021) 52.
  * Rincon [2019] F. Rincon,  Dynamo theories,  J. Plasma Phys 85 (2019).
  * Dormy et al. [2004] E. Dormy, A. M. Soward, C. A. Jones, D. Jault, P. Cardin,  The onset of thermal convection in rotating spherical shells,  J. Fluid Mech. 501 (2004) 43–70.
  * Iyer et al. [2020] K. Iyer, J. Scheel, J. Schumacher, K. Sreenivasan,  Classical 1/3 scaling of convection holds up to $\text{Ra}=10^{15}$,  Proc. Natl. Acad. Sci. USA 117 (2020) 7594–7598.
  * Wicht [2002] J. Wicht, Inner-core conductivity in numerical dynamo simulations, 2002.
  * Forte et al. [1995] A. M. Forte, J. X. Mitrovica, R. L. Woodward,  Seismic-geodynamic determination of the origin of excess ellipticity of the core-mantle boundary,  Geophys. Res. Lett. 22 (1995) 1013–1016.
  * Yoshida et al. [1996] S. Yoshida, I. Sumita, M. Kumazawa,  Growth model of the inner core coupled with the outer core dynamics and the resulting elastic anisotropy,  J. Geophys. Res. 101 (1996) 28085–28103.
  * Chongsiripinyo [2019] K. Chongsiripinyo, Decay of stratified turbulent wakes behind a bluff body, University of California, San Diego, 2019.
  * Taylor [2008] J. Taylor, Numerical simulations of the stratified oceanic bottom boundary layer publication date, 2008\. URL: https://escholarship.org/uc/item/5s30n2ts.
  * Brown et al. [2000] P. Brown, R. Falgout, J. Jones, S. Comput, Semicoarsening multigrid on distributed memory machines, 2000\. URL: http://www.siam.org/journals/sisc/21-5/33914.html.
  * Falgout and Jones [2000] R. Falgout, J. Jones, Multigrid on massively parallel architectures, 2000\. URL: http://www.llnl.gov/tid/Library.html.
  * Falgout et al. [2002] R. Falgout, J. Jones, U. Yang, The design and implementation of hypre, a library of parallel high performance preconditioners, 2002\.
  * Ruge and Stüben [1987] J. Ruge, K. Stüben, Algebraic multigrid, 1987. URL: http://www.siam.org/journals/ojsa.php.
  * Long et al. [2020] R. Long, J. Mound, C. J. Davies, S. M. Tobias,  Thermal boundary layer structure in convection with and without rotation,  Phys. Rev. Fluids. 5 (2020).
  * Tennekes and Lumley [1972] H. Tennekes, J. Lumley, A First Course in Turbulence, The MIT press, 1972\.
  * Brucker and Sarkar [2010] K. Brucker, S. Sarkar,  A comparative study of self-propelled and towed wakes in a stratified fluid,  J. Fluid Mech. 652 (2010) 373–404.
  * Kooij et al. [2018] G. L. Kooij, M. A. Botchev, E. M. Frederix, B. J. Geurts, S. Horn, D. Lohse, E. P. van der Poel, O. Shishkina, R. Stevens, R. Verzicco,  Comparison of computational codes for direct numerical simulations of turbulent rayleigh–bénard convection,  Computers and Fluids 166 (2018) 1–8.
  * Busse et al. [1998] F. Busse, E. Grote, A. Tilgner,  On convection driven dynamos in rotating spherical shells,  Studia Geophysica et Geodaetica 42 (1998) 211.
  * Christensen et al. [1998] U. Christensen, P. Olson, G. A. Glatzmaier,  A dynamo model interpretation of geomagnetic field structures,  Geophysical Research Letters 25 (1998) 1565–1568.
  * Christensen et al. [1999] U. Christensen, P. Olson, G. Glatzmaier,  Numerical modelling of the geodynamo: a systematic parameter study,  Geophysical Journal International 138 (1999) 393–409.
  * Dormy et al. [1998] E. Dormy, P. Cardin, D. Jault,  Mhd flow in a slightly differentially rotating spherical shell, with conducting inner core, in a dipolar magnetic field,  Earth and Planetary Science Letters 160 (1998) 15–30.
  * Glatzmaier [1984] G. A. Glatzmaier,  Numerical simulations of stellar convective dynamos. i. the model and method,  Journal of Computational Physics 55 (1984) 461–484.
  * Sakuraba and Kono [1999] A. Sakuraba, M. Kono,  Effect of the inner core on the numerical solution of the magnetohydrodynamic dynamo,  Physics of the Earth and Planetary Interiors 111 (1999) 105–121.
  * Tilgner [1999] A. Tilgner,  Spectral methods for the simulation of incompressible flows in spherical shells,  International journal for numerical methods in fluids 30 (1999) 713–724.
  * Kageyama et al. [1995] A. Kageyama, T. Sato, C. S. Groupa),  Computer simulation of a magnetohydrodynamic dynamo. ii,  Physics of Plasmas 2 (1995) 1421–1431.