A direction splitting scheme for Navier-Stokes-Boussinesq system in spherical shell geometries

TL;DR

This paper presents a second-order direction-splitting numerical scheme for solving the Navier-Stokes-Boussinesq system in spherical shells, utilizing overset Yin-Yang grids and parallel domain decomposition for improved stability and scalability.

Contribution

It introduces a novel second-order direction-splitting method on overset Yin-Yang grids for spherical shell geometries, avoiding pole singularities and enabling parallel computation.

Findings

Method achieves second-order temporal accuracy.

Scheme demonstrates stability and weak scalability.

Validated on manufactured and Landau solutions.

Abstract

This paper introduces a formally second-order direction-splitting method for solving the incompressible Navier-Stokes-Boussinesq system in a spherical shell region. The equations are solved on overset Yin-Yang grids, combined with spherical coordinate transforms. This approach allows to avoid the singularities at the poles and keeps the grid size relatively uniform. The downside is that the spherical shell is subdivided into two equally sized, overlapping subdomains that requires the use of Schwarz-type iterations. The temporal second order accuracy is achieved via an Artificial Compressibility (AC) scheme with bootstrapping. The spatial discretization is based on second order finite differences on the Marker-And-Cell (MAC) stencil. The entire scheme is implemented in parallel using a domain decomposition iteration, and a direction splitting approach for the local solves. The stability, accuracy and weak scalability of the method is verified on a manufactured solution of the Navier-Stokes-Boussinesq system and on the Landau solution of the Navier-Stokes equations on the sphere.

Equations

$$\frac{\partial \mathbf{u}}{\partial t}+(\mathbf{u}\cdot \nabla )\mathbf{u}+\nabla p-\mathrm{Pr}\Delta \mathbf{u}$$

$$\nabla \cdot \mathbf{u}$$

$$\mathbf{u}$$

$$\frac{\partial T}{\partial t}+(\mathbf{u}\cdot \nabla )T-\Delta T$$

$$T$$

$$\Omega =\left\{\left(r,\theta ,\varphi \right)\in [R_1,R_2]\times \left[0,\pi \right]\times \left[0,2\pi \right)\right\}.$$

$$\{\begin{array}{ll}x& =r\sin \theta \cos \varphi \\ y& =r\sin \theta \sin \varphi \\ z& =r\cos \theta ,\end{array}$$

$${\Omega }_1:=\left\{\left(r,\theta ,\varphi \right)\in [R_1,R_2]\times \left[\frac{\pi }{4}-\epsilon ,\frac{7\pi }{4}+\epsilon \right]\times \left[\frac{\pi }{4}-\epsilon ,\frac{3\pi }{4}+\epsilon \right]\right\},$$

$$\{\begin{array}{ll}x& =-r\sin \tilde{\theta }\cos \tilde{\varphi }\\ y& =r\cos \tilde{\theta }\\ z& =r\sin \tilde{\theta }\sin \tilde{\varphi },\end{array}$$

$${\Omega }_2:=\left\{\left(r,\tilde{\theta },\tilde{\varphi }\right)\in [R_1,R_2]\times \left[\frac{\pi }{4}-\epsilon ,\frac{7\pi }{4}+\epsilon \right]\times \left[\frac{\pi }{4}-\epsilon ,\frac{3\pi }{4}+\epsilon \right]\right\}.$$

$${\partial }_{t}T-\kappa \Delta T$$

$$T$$

$$\Delta ={\mathrm{D}}_{rr}+{\mathrm{D}}_{\theta \theta }+{\mathrm{D}}_{\varphi \varphi },{\mathrm{D}}_{rr}:=\frac{1}{r^2}{\partial }_r\left(r^2{\partial }_r\right),{\mathrm{D}}_{\theta \theta }:=\frac{1}{r^2\sin \theta }{\partial }_{\theta }\left(\sin \theta {\partial }_{\theta }\right),\text{ and }{\mathrm{D}}_{\varphi \varphi }:=\frac{{\partial }_{\varphi \varphi }}{r^2{\sin }^2\theta }.$$

$$\left[\mathrm{I}-\frac{\tau }{2}{\mathrm{D}}_{rr}\right]\left[\mathrm{I}-\frac{\tau }{2}{\mathrm{D}}_{\theta \theta }\right]\left[\mathrm{I}-\frac{\tau }{2}{\mathrm{D}}_{\varphi \varphi }\right]\frac{\delta T^{n+1}}{\tau }=\Delta T^n,$$

$${\hat{\mathrm{D}}}_{\theta \theta }:=\frac{1}{R_1^2\sin \theta }{\partial }_{\theta }\left(\sin \theta {\partial }_{\theta }\right),{\hat{\mathrm{D}}}_{\varphi \varphi }:=\frac{{\partial }_{\varphi \varphi }}{R_1^2{\sin }^2{\theta }_1},\text{ and }\hat{\Delta }:={\mathrm{D}}_{rr}+{\hat{\mathrm{D}}}_{\theta \theta }+{\hat{\mathrm{D}}}_{\varphi \varphi },$$

$$-{\left({\hat{\mathrm{D}}}_{\theta \theta }T,T\right)}_{\omega }\ge 0,-{\left({\hat{\mathrm{D}}}_{\varphi \varphi }T,T\right)}_{\omega }\ge 0$$

$$-{\left(\left[{\hat{\mathrm{D}}}_{\theta \theta }-{\mathrm{D}}_{\theta \theta }\right]T,T\right)}_{\omega }\ge 0\text{ and }-{\left(\left[{\hat{\mathrm{D}}}_{\varphi \varphi }-{\mathrm{D}}_{\varphi \varphi }\right]T,T\right)}_{\omega }\ge 0.$$

$${\left(-\hat{\Delta }T,T\right)}_{\omega }\ge {\left(-\Delta T,T\right)}_{\omega }.$$

$$\frac{\delta T^{n+1}}{\tau }=\Delta T^{\ast ,n+1/2},$$

$$\left[\mathrm{I}-\frac{\tau }{2}{\mathrm{D}}_{rr}\right]\left[\mathrm{I}-\frac{\tau }{2}{\hat{\mathrm{D}}}_{\theta \theta }\right]\left[\mathrm{I}-\frac{\tau }{2}{\hat{\mathrm{D}}}_{\varphi \varphi }\right]\frac{\delta T^{n+1}}{\tau }=\Delta T^{\ast ,n+1/2}-\frac{1}{2}\hat{\Delta }\delta T^n.$$

$$\frac{{\tau }^2}{2}\hat{\Delta }\frac{{\delta }^2{T}^{n+1}}{{\tau }^2}+\left[\frac{{\tau }^2}{4}({\mathrm{D}}_{rr}{\hat{\mathrm{D}}}_{\theta \theta }+{\mathrm{D}}_{rr}{\hat{\mathrm{D}}}_{\varphi \varphi }+{\hat{\mathrm{D}}}_{\theta \theta }{\hat{\mathrm{D}}}_{\varphi \varphi })-\frac{{\tau }^3}{8}{\mathrm{D}}_{rr}{\hat{\mathrm{D}}}_{\theta \theta }{\hat{\mathrm{D}}}_{\varphi \varphi }\right]\frac{\delta T^{n+1}}{\tau }.$$

$$\tau \sum _{n=1}^{N-1}\frac{\Vert T^{n+1}-T^n{\Vert }_{\omega }^2}{{\tau }^2}$$

$$\le \frac{1}{2}\Vert \nabla T^1{\Vert }_{\omega }^2+\frac{1}{4}\left(\Vert {\partial }_{\theta }\left(T^1-T^0\right){\Vert }_{{\omega }_1}^2+\Vert {\partial }_{\varphi }\left(T^1-T^0\right){\Vert }_{{\omega }_2}^2\right),$$

$$[\mathrm{I}-\frac{\tau }{2}\hat{\Delta }$$

$$\frac{\delta T^{n+1}}{\tau }$$

$$+\left[\frac{{\tau }^2}{4}\left({\mathrm{D}}_{rr}{\hat{\mathrm{D}}}_{\theta \theta }+{\mathrm{D}}_{rr}{\hat{\mathrm{D}}}_{\varphi \varphi }+{\hat{\mathrm{D}}}_{\theta \theta }{\hat{\mathrm{D}}}_{\varphi \varphi }\right)-\frac{{\tau }^3}{8}{\mathrm{D}}_{rr}{\hat{\mathrm{D}}}_{\theta \theta }{\hat{\mathrm{D}}}_{\varphi \varphi }\right]\frac{\delta T^{n+1}}{\tau }=0.$$

$$-\frac{1}{2}{\left(\left[\Delta -\hat{\Delta }\right]\left(T^{n+1}-2T^n+T^{n-1}\right),T^{n+1}-T^n\right)}_{\omega }$$

$$=\frac{1}{4}\left[\Vert {\partial }_{\theta }\left(T^{n+1}-T^n\right){\Vert }_{{\omega }_1}^2-\Vert {\partial }_{\theta }\left(T^n-T^{n-1}\right){\Vert }_{{\omega }_1}^2+\Vert {\partial }_{\theta }\left(T^{n+1}-2T^n+T^{n-1}\right){\Vert }_{{\omega }_1}^2\right]$$

$$+\frac{1}{4}\left[\Vert {\partial }_{\varphi }\left(T^{n+1}-T^n\right){\Vert }_{{\omega }_2}^2-\Vert {\partial }_{\varphi }\left(T^n-T^{n-1}\right){\Vert }_{{\omega }_2}^2+\Vert {\partial }_{\varphi }\left(T^{n+1}-2T^n+T^{n-1}\right){\Vert }_{{\omega }_2}^2\right].$$

$$-{\left(\Delta T^{n+1/2},T^{n+1}-T^n\right)}_{\omega }=\frac{1}{2}\left(\Vert \nabla T^{n+1}{\Vert }_{\omega }^2-\Vert \nabla T^n{\Vert }_{\omega }^2\right).$$

$${\left({\mathrm{D}}_{rr}{\hat{\mathrm{D}}}_{\theta \theta }\delta T^{n+1},\delta T^{n+1}\right)}_{\omega }^2=\underset{\Omega }{\int }\frac{r^2\sin \theta }{R_1^2}|{\partial }_{r\theta }\delta T^{n+1}{|}^2,$$

$${\left({\mathrm{D}}_{rr}{\hat{\mathrm{D}}}_{\varphi \varphi }\delta T^{n+1},\delta T^{n+1}\right)}_{\omega }=\underset{\Omega }{\int }\frac{r^2\sin \theta }{R_1^2{\sin }^2{\theta }_1}|{\partial }_{r\varphi }\delta T^{n+1}{|}^2,$$

$${\left({\hat{\mathrm{D}}}_{\theta \theta }{\hat{\mathrm{D}}}_{\varphi \varphi }\delta T^{n+1},\delta T^{n+1}\right)}_{\omega }=\underset{\Omega }{\int }\frac{r^2\sin \theta }{R_1^4{\sin }^2{\theta }_1}|{\partial }_{\theta \varphi }\delta T^{n+1}{|}^2,$$

$$-{\left({\mathrm{D}}_{rr}{\hat{\mathrm{D}}}_{\theta \theta }{\hat{\mathrm{D}}}_{\varphi \varphi }\delta T^{n+1},\delta T^{n+1}\right)}_{\omega }=\underset{\Omega }{\int }\frac{r^2\sin \theta }{R_1^4{\sin }^2{\theta }_1}|{\partial }_{r\theta \varphi }\delta T^{n+1}{|}^2.$$

Full text

††journal: ArXiv.org

1 Introduction

This article presents a new direction-splitting scheme for solving the incompressible Navier-Stokes-Boussinesq system:

[TABLE]

[TABLE]

in a spherical shell domain that can be defined in terms of a spherical coordinate triple $\left(r,\theta,\phi\right)$ as:

[TABLE]

In the above, ${\bf g}$ is the unit vector in the direction of gravity, and $\Pr,\mathrm{Ra}$ are the Prandtl and Rayleigh numbers, respectively. The system (1.1)-(1.2) models the flow of a heat conducting fluid, under the assumption that the temperature-induced density variation influences significantly only the buoyancy force and the fluid remains incompressible. It is widely applied to model the flow in the atmospheric boundary layer ([3]), oceanic flows ([4]), as well as, if combined with an equation for the magnetic field, the flow in the Earth’s dynamo ([5]). Even though for the most part of the discussion, we assume homogeneous Dirichlet boundary conditions on the two spherical surfaces $r=R_{1},r=R_{2}$, the approach is applicable to Neumann and Robin boundary conditions as well.

One widely used approach for numerical approximations of differential equations in spherical shell geometries is based on the use of a spherical transformation that transforms the domain into a parallelepiped. The obvious advantage of this approach is the simple computational domain, which allows for the use of structured grids and the efficient schemes developed for them. Moreover, the grid can naturally follow the geometry of the domain, without requiring too many cells, as would possibly be in the case of a Cartesian formulation. However, the singularity of the transformation and the grid convergence near the poles have for many years been a difficulty in the development of accurate finite difference and pseudo-spectral schemes. Several different treatments have been proposed for dealing with these problems. For example, in [6], the pole singularity issue is avoided by replacing the equations at the poles with equations, analogous to boundary conditions, while in [7], a redefinition of the singular coordinates is proposed. Other suggested approaches include applying L’Hospital’s rule [8] to singular terms and switching to Cartesian formulations around the poles [9]. On the other hand, the grid convergence has been a more serious problem. In particular, it produces a solution with uneven resolution, requires very small time steps for explicit or IMEX schemes since the time step size is limited by the minimum grid size, and causes convergence problems for iterative solvers. Therefore, different grid systems have been suggested in the literature that give quasi-uniform resolution and avoid the grid convergence problems. One such approach is the ”cubed sphere” of [10], which is a grid that covers a spherical surface with six components corresponding to six faces of a cube. Even though, the resulting grid is quasi-uniform, it still has singularities at the corner points of the faces and it is non-orthogonal. Some of the other suggested unstructured grids include the isocahedral grid of [11] and non-orthogonal rhombahedral grid of [12].

In this study we adopt an alternative approach, proposed by [13], employing the so-called Yin-Yang grids. It starts with a decomposition of the domain into two overlapping subdomains, combined with two different spherical transforms whose axes are perpendicular to each other, cf. Fig. 1. As a result, both subdomains are transformed into identical parallelepipeds that can be gridded with the same uniform grids. This approach automatically removes the transforms singularities at the poles, at the expense of the introduction of two subdomains, so that the two local solutions must be coupled by means of Schwarz-type iterations. It has been used for simulations of mantle convection [14], core collapse supernovae [15], atmospherical general circulation model [16] and visualization in spherical regions [17]. Some advantages of Yin-Yang approach are that the metric tensors are simple, the resolution is quasi-uniform, and it requires modest programming effort for extending the code from a single latitude-longitude grid. The main novelty of this paper is that the Yin-Yang domain decomposition is combined with a direction splitting time discretization that, in case of linear parabolic equations, is unconditionally stable on grids on the spherically transformed domains. The advection can be included either in an IMEX fashion, or by including the linearized advection operator into the entire operator that is further split direction-wise. The resulting splitting scheme is conditionally stable, since the direction-wise operators are not positive, but our numerical experience demonstrated that the second approach yields an algorithm that has better stability performance. This is why the rest of the paper concerns only this type of schemes. To our knowledge, the stability of the direction splitting approach has not been rigorously studied in the context a spherical coordinate system. Therefore, we prove below that it is unconditionally stable in case of a scalar heat equation, in a simply shaped domain (in terms of spherical coordinates). The case of the full Navier-Stokes-Boussinesq system is more involved and we do not provide a rigorous proof here. However, our numerical experience shows that the a direction splitting is still unconditionally stable if the advection terms are omitted and if the velocity-pressure decoupling is done via the AC method proposed in [1].

The rest of the paper is organized as follows. In the next section, we briefly recall the definition of the Yin-Yang domain decompostion. In Section $3$ we present the numerical scheme for the advection-diffusion and Navier-Stokes equations on each of the subdomains. In Section $4$, we discuss the implementation details, and in Section $5$ we present the numerical experiments.

2 Spatial discretization and the Yin-Yang grid

In this section, we briefly recall the definition of the composite Yin-Yang grid following [13]. The grid consists of two identical overlapping latitude-longitude grids whose axes are perpendicular to each other. The Yin grid is based on a spherical transformation

[TABLE]

and covers the region

[TABLE]

where $\varepsilon\ll 1$ is a parameter determining the overlap. The Yang grid is obtained via another spherical transformation:

[TABLE]

such that its axes is perpendicular to the axes of the Yin transform, and covers the region

[TABLE]

The choice of the second axes should be such that the Yang grid fully covers the gap of the Yin one, and the overlapping subregions are of the same size (see Fig. 1,). Otherwise, it is identical to the Yin grid modulo two rotations. The resulting Yin-Yang grids are quasiuniform, the coordinate transformations from $(r,\theta,\phi)$ to $(r,\tilde{\theta},\tilde{\phi})$ and its inverse, as well as the metric tensors on both grids are identical. As a consequence, the methods and codes developed for the standard latitude-longitude grid can be applied to both grids.

3 Direction-splitting time discretization

3.1 Preliminaries

In the sequel of the paper we will frequently make use of the following notations. For a time sequence $w^{k},k=1,2,\dots$ we denote the average between two time levels as $w^{k+1/2}=(w^{k+1}+w^{k})/2$, and the explicit extrapolation to level $k+1/2$ by $w^{*,k+1/2}=3w^{k}/2-w^{k-1}/2$. For two regular enough functions $u,v$ defined in the spherical shell we denote their weighted $L^{2}$ inner product as: $(u,v)_{\omega}:=\int\limits_{R_{1}}^{R_{2}}\int\limits_{0}^{\pi}\int\limits_{0}^{2\pi}uv\omega drd\theta d\phi,$ where $\omega$ denotes a non-negative weight. In most cases the weight is given by the weight of the spherical transform $\omega=r^{2}\sin\theta$, however, in some of the estimates given below, the weight will be appropriately modified. The corresponding norm is given by $\|u\|_{\omega}^{2}=(u,u)_{\omega}$.

3.2 Direction splitting of the advection-diffusion equation

Since the PDEs are identical in both domains, it is sufficient to develop the numerical scheme for the Yin domain. Then the Schwarz domain decomposition method can be used to iterate between the subdomains.

We first present Douglas [18] type direction splitting scheme for the heat equation. Consider

[TABLE]

where the Laplacian in spherical coordinates is given by

[TABLE]

The Douglas direction splitting scheme for this equation can be summarized in the following factorized form:

[TABLE]

where $\delta T^{n+1}:=T^{n+1}-T^{n}$ denotes the first time difference of the time sequence $T^{k}$, $\tau$ is the time step, and $\mathrm{I}$ is the identity operator. We first notice that this splitting can be considered as an Euler explicit scheme whose time difference operator is multiplied by $\left[\mathrm{I}-\frac{\tau}{2}\mathrm{D}_{rr}\right]\left[\mathrm{I}-\frac{\tau}{2}\mathrm{D}_{\theta\theta}\right]\left[\mathrm{I}-\frac{\tau}{2}\mathrm{D}_{\phi\phi}\right]$, that is a consistent perturbation of $\mathrm{I}$ and stabilizes the scheme. If the spatial derivative operators are positive and commute with respect to some inner product, the stability of this scheme is not hard to establish.

Unfortunately, $\mathrm{D}_{rr},\mathrm{D}_{\theta\theta},$ and $\mathrm{D}_{\phi\phi}$ do not commute with respect to the weighted product $(.,.)_{\omega}$, and their positivity is far from being clear. The main obstacle to commutativity of the one-dimensional operators comes from the non-constant terms in the denominators of $\mathrm{D}_{\theta\theta}$ and $\mathrm{D}_{\phi\phi}$. Therefore, the scheme (3.2) should be modified as follows. We first introduce the modified spatial operators:

[TABLE]

where $\displaystyle\theta_{1}=\frac{\pi}{4}-\varepsilon$. Then it is easy to check that $\mathrm{D}_{rr},\hat{\mathrm{D}}_{\theta\theta}$ and $\hat{\mathrm{D}}_{\phi\phi}$, supplied with zero Dirichlet boundary conditions, commute. Moreover,

[TABLE]

and

[TABLE]

These inequalities immediately yield that:

[TABLE]

In order to obtain an unconditionally stable second order scheme, we start from the second order Adams-Bashforth scheme:

[TABLE]

and stabilize it by multiplying the time difference in the left hand side by $\left[\mathrm{I}-\frac{\tau}{2}\mathrm{D}_{rr}\right]\left[\mathrm{I}-\frac{\tau}{2}\hat{\mathrm{D}}_{\theta\theta}\right]\left[\mathrm{I}-\frac{\tau}{2}\hat{\mathrm{D}}_{\phi\phi}\right].$ Since this perturbation is only first order consistent with the identity operator, we subtract from the right hand side the first order perturbation term, taken at the previous time level. The resulting splitting scheme reads:

[TABLE]

Note that, assuming enough regularity of the exact solution in space and time, this is a second-order perturbation of the second-order explicit Adams-Bashforth scheme (3.2), the perturbation being given by:

[TABLE]

We have the following stability result for the scheme (3.6).

Theorem 3.1**.**

Assuming enough regularity of the exact solution $T$ of the semi-discrete scheme (3.6), it is unconditionally stable; more precisely, it satisfies the following estimate:

[TABLE]

where $\omega_{1}=\left(1-\frac{r^{2}}{R^{2}_{1}}\right)\sin\theta\geq 0$ and $\omega_{2}=\left(\frac{r^{2}}{R^{2}_{1}}-1\right)\frac{\sin\theta}{\sin^{2}\theta_{1}}\geq 0$.

Proof.

Expanding the left hand side of (3.6) we get:

[TABLE]

Rearranging all the $\Delta$ and $\hat{\Delta}$ terms, we obtain

[TABLE]

Next we multiply (3.9) by $v=\delta T^{n+1}$ and integrate by parts. Then the second term gives

[TABLE]

The third term is

[TABLE]

The remaining terms are all dissipative:

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

Substituting (3.10)-(3.15) into (3.9), and summing for $n=1,\dots,N-1$ completes the proof. ∎

The factorized scheme for the advection-diffusion equation (1.2) is obtained in a similar fashion and takes the following form:

[TABLE]

3.3 Direction-splitting discretization of the Navier-Stokes system

Now we present the direction splitting scheme for the Navier-Stokes equations (1.1). Our numerical scheme is based on the AC regularization:

[TABLE]

where $\chi=\mathcal{O}\left(1\right)$ is an artificial compressibility regularization parameter, and $\mathrm{Re}$ is the Reynolds number. It is well-known that the resulting approximation $\left(\mathbf{u}_{1},p_{1}\right)$ is first-order accurate in time (see [19]). A second order scheme can be constructed using the bootstrapping approach of [1, 2], which requires additionally to solve the system:

[TABLE]

$p_{1}$ being given by (3.17). In the following, for the sake of brevity, we will only discuss the direction splitting implementation of the first order approximation (3.17). The higher order correction for $\mathbf{u}_{2},p_{2}$ is solved identically. First, consider the standard semi-implicit Crank-Nicholson approximation of the system for $\left(\mathbf{u}_{1},p_{1}\right)$:

[TABLE]

Note that we use the second order velocity $\mathbf{u}_{2}$ as advecting velocity, which allows us to assemble a single linear system for both systems. We can rewrite the momentum equation by eliminating $p^{n+1}_{1}$ from the first equation:

[TABLE]

In order to produce a factorized scheme for each velocity component, the $\nabla\nabla\!{\cdot}$ operator must be also split somehow, and we use the Gauss-Seidel type splitting of the $\nabla\nabla\!{\cdot}$ operator, which was originally proposed in [2] in the Cartesian case:

[TABLE]

3.3.1 Equation for the $r$-component of the velocity

Using the mass conservation equation $\nabla\cdot\mathbf{u}=0$, it is possible to write the first component of the system as follows:

[TABLE]

where $\mathbf{u}\cdot\nabla v=u_{r}\mathrm{\partial}_{r}v+u_{\theta}\frac{\mathrm{\partial}_{\theta}v}{r}+u_{\phi}\frac{\mathrm{\partial}_{\phi}v}{r\sin\theta}$ is the advection operator. Let $\mathrm{L}_{rr},\mathrm{L}_{r\theta}$ and $\mathrm{L}_{r\phi}$ be the differential operators that act in each space direction:

[TABLE]

[TABLE]

The factorized scheme for the $r$-component takes the following form:

[TABLE]

3.3.2 Equation for the $\theta$–component of the velocity

Again using $\nabla\cdot\mathbf{u}=0$, the $\theta$-component of the momentum equation can be expressed as:

[TABLE]

Let $\mathrm{L}_{\theta r},\mathrm{L}_{\theta\theta}$ and $\mathrm{L}_{\theta\phi}$ be defined as follows:

[TABLE]

[TABLE]

[TABLE]

The factorized scheme for the $\theta$-component takes the following form:

[TABLE]

3.3.3 Equation for the $\phi$–component of the velocity

The $\phi$-component of the momentum equation is given by:

[TABLE]

Let $\mathrm{L}_{\phi r},\mathrm{L}_{\phi\theta}$ and $\mathrm{L}_{\phi\phi}$ be defined as follows:

[TABLE]

[TABLE]

The factorized scheme for the $\phi$-component is then:

[TABLE]

3.3.4 Pressure update

[TABLE]

4 Implementation and parallelization

The equations (3.16), (3.19)-(3.21) are solved as a sequence of 1D equations in each space direction. For example, solving (3.16) consists of the following steps:

[TABLE]

Similar strategy is applied for the Navier-Stokes approximation. Each 1D system is spatially approximated using second-order centered finite differences on a non-uniform grid. In order to ensure the inf-sup stability, the unknowns are approximated on a MAC grid, where the velocity components are stored at the face centers of the cells, while the scalar variables are stored at the cell centers.

To solve the system on each domain in parallel we use the approach developed in [20], where we first perform Cartesian domain decomposition of both computational grids using MPI, and then solve the resulting set of tridiagonal linear systems using domain-decomposition-induced Schur complement technique. Note, that the Schur complement can be computed explicitly (see [20] for details) and so the system in each direction can be solved directly by the Thomas algorithm, avoiding the need of iterations on each of the two subdomain. Then, in order to obtain the approximation on the entire spherical shell, we iterate between the Yin and Yang grids using either additive or multiplicative overlapping Schwarz methods. The solution on each grid is computed using only boundary data that is interpolated from the currently available solution on the other grid, using Lagrange interpolation.

In the additive Schwarz implementation, we use an even total number of CPUs. Then we split the global communicator into two equal parts, and assign to each grid one of the communicators. In the multiplicative Schwarz implementation, we use the global communicator to solve the problem on each grid sequentially.

The overall solution procedure in case of the multiplicative Schwarz iteration can be summarized as follows:

Algorithm 4.1**.**

Repeat until convergence:

• For $i=1,2$

Obtain interpolated boundary values $T_{bd}$ for $\partial\Omega_{i}$ from $\Omega_{3-i}$.

• 2)

Solve the temperature equation in $\Omega_{i}$ with using extrapolated velocity values $\mathbf{u}_{2}^{*,n+1/2}$.

• 3)

Obtain interpolated boundary values $\mathbf{u}_{bd}$ for $\partial\Omega_{i}$ from $\Omega_{3-i}$.

• 4)

If \left|\bigintsss\limits_{\partial\Omega_{i}}\mathbf{u}_{bd}\cdot\mathbf{n}\right|\geq\mathrm{tol}, then minimize the functional ($\varepsilon\ll 1$):

[TABLE]

using the Conjugate Gradient Algorithm until $\mathrm{J}\left(\cdot\right)\leq\mathrm{tol}$.

• 5)

Update $\mathbf{u}_{bd}:={\bf v}$ and solve the momentum equation in $\Omega_{i}$ with the interpolated Dirichlet boundary conditions in $\theta,\phi$ directions and with the original boundary conditions in the $r$ direction.

• 6)

Compute the pressure in $\Omega_{i}$ using the second equation in (3.17).

• 7)

Interpolate the pressure values at the boundary of $\partial\Omega_{i}$ using the available pressure on $\Omega_{3-i}$.

• End for.

Step 4 is meant to ensure that there is no spurious mass flux generated through the internal (artificial) boundaries due to the interpolation. It is optional, and as our numerical experience shows, it rarely changes significantly the results. Therefore, it is skipped while producing the numerical results presented in the next section. Skipping Step 7, however, can seriously reduce the rate of convergence of the Schwarz iteration, as observed in the numerical simulations. Clearly, the AC method for the Navier-Stokes equations does not require boundary conditions on the pressure. Nevertheless, the exchange of the pressure values does influence the pressure gradient that appears in (3.19)-(3.21), and thus it seems to influence significantly the convergence of the overall iteration. This effect is not well understood and while some other authors (see for example [21]) also interpolate the pressure values near the internal boundaries, others (e.g. [22] ) interpolate only the velocity on the internal boundaries.

Another interesting feature of the domain decomposition iteration described above is that it allows to use the previously computed iterates in order to reduce the splitting error of the direction splitting approximation. For example, if the factorized form of the direction-split approximation for a given quantity $\psi$ is given by:

[TABLE]

where the superscript $n$ denotes the time level of the solution and $k$ denotes the domain decomposition iteration level, then the splitting error can be reduced by using the modified equation:

[TABLE]

where $L_{\psi}=L_{\psi,r}+L_{\psi,\theta}+L_{\psi,\phi}$. Indeed, in (4.1), the splitting error term $(L_{\psi,r}L_{\psi,\theta}+L_{\psi,r}L_{\psi,\phi}+L_{\psi,\theta}L_{\psi,\phi}-L_{\psi,r}L_{\psi,\theta}L_{\psi,\phi})(\psi^{n,k}-\psi^{n-1})$ at iteration level $k$ has been reduced by the same term at the previous iteration level $(L_{\psi,r}L_{\psi,\theta}+L_{\psi,r}L_{\psi,\phi}+L_{\psi,\theta}L_{\psi,\phi}-L_{\psi,r}L_{\psi,\theta}L_{\psi,\phi})(\psi^{n,k-1}-\psi^{n-1})$. If this error reduction is employed, then the iteration becomes a block-preconditioned overlapping domain decomposition iteration, the preconditioner being the factorized operator $(I-L_{\psi,r})(I-L_{\psi,\theta})(I-L_{\psi,\phi})$. We must also remark here that this iteration needs to converge to an accuracy of the order of $\tau^{2}$ for the solution of equation (3.17) and $\tau^{3}$ for the solution of equation (3.18), in order to preserve the second order accuracy of the overall algorithm.

5 Numerical tests

5.1 Time and space convergence

We verify the convergence rates in space and time using the following manufactured solution, given in a Cartesian form:

[TABLE]

The parameters used in this test are $R_{1}=1,R_{2}=2,Ra=1,Pr=1$, and the grids used in the tests are uniform in each direction. The convergence of the approximation is tested using both, the additive and multiplicative versions of the scheme. The grid used for the time convergence tests consists of $20\times 92\times 192$ MAC cells on each of the two subdomains. The solution error is computed at the final time $T_{f}=10$. For the space convergence tests, the time step is chosen small enough to not influence the overall error, $\tau=0.0001$, and the final time is $T_{f}=1$. The grid diameter is computed as the maximum diameter of the MAC cells in Cartesian coordinates. In both cases the domain decomposition iterations are converged so that the $l^{2}$ norm of the difference between two subsequent iterates, for any of the computed quantities is less than $10^{-6}$ ($l^{2}$ norm denotes the standard mid-point approximation to the $L^{2}$ norm). Also, the splitting error reduction, as outlined by equation (4.1), is employed at each iteration.

The graphs of the $l^{2}$ norm of the errors in both cases are presented in Figure 2. They clearly demonstrate the second order accuracy of the scheme in space and time.

Next we verify the accuracy of the proposed algorithm on a physically more relevant analytic solution of the Navier-Stokes equations in a spherical setting, due to Landau (see [23] and [24] for a recent review). The source term of the equations is equal to zero in this case, and the solution is steady and axisymmetric. In all cases presented in figure 4 the multiplicative Schwarz version of the algorithm is used with its convergence tolerance being set to $10^{-6}$, the time step is equal to $10^{-3}$, and $R_{1}=1,R_{2}=2$. We first present in the top left graph of figure 4 the $l^{2}$ error for the velocity and pressure as a function of the grid diameter, at $Re=1$ and the overlap is $\epsilon=0.1$ The scheme clearly exhibits again a second order convergence rate in space. In the top right graph we demonstrate the influence of the overlap size on the error at $Re=1$, the grid size in the $r,\phi,\theta$ directions being $2.7778\times 10^{-2},1.7027\times 10^{-2},3.6121\times 10^{-2}$ correspondingly. The effect of the overlap on the error is insignificant, however, it seriously impacts the stability of the algorithm i.e. the increase of the overlap improves the stability, particularly at large Reynolds numbers. Finally, the bottom graph demonstrates the effect of the Reynolds number on the error. Again, the overlap is $0.1$ and the grid sizes are equal to $2.7778\times 10^{-2},1.7027\times 10^{-2},3.6121\times 10^{-2}$ . We should note that the exact solution for the velocity scales like $Re^{-1}$ and therefore the errors in the graph are multiplied by the corresponding Reynolds number. Clearly, the oscillations in the error decrease slower with the increase of the Reynolds number. These oscillations are due to the artificial compressibility algorithm, since the initial data for the pressure corresponds to a divergence-free velocity, while the pressure evolution is determined by a perturbed continuity equation (see [25] and [26] for a detailed discussion on this issue).

5.2 Weak parallel efficiency

Next we test the parallel efficiency of the code based on the scheme introduced in the previous section. Since we are interested in solving large size problems, we only measure the weak scalability of our code. The problem size is $100\times 100\times 100$ grid cells per each of the Yin and Yang grids on each CPU, and the maximum number of CPUs used is $960$. Besides, since in the possible applications of this technique (atmospheric boundary layer, Earth’s dynamo) the thickness of the spherical shell is much smaller than the diameter of the shell, we use a two-dimensional grid of processors for the grid partitioning. It must be noted though, that making the grid partitioning three dimensional does not change much the parallel efficiency results presented in this section. The scaling efficiency is computed as the ratio of the CPU time on 32 cores divided by the CPU time on $n\geq 32$ cores. The reason for this definition of efficiency is that the particular cluster used in the scaling tests has processors containing 32 cores each, and the efficiency drops very significantly between 1 and 32 cores (to about 75%). After this, when the number of cores is a multiple of 32 the efficiency remains very close to the one at 32 cores. One possible explanation of this phenomenon is that in case when the number of cores is significantly less than 32 cores, they need to share the memory bandwidth and cache with a smaller number of cores, since presumably the rest of the available cores on the given processor are idle (see e.g. [27], p. 152). Again, we are interested in very large computations, and therefore, using a minimum of 32 cores is very reasonable.

The scaling results are performed using the Compute Canada (see https://www.computecanada.ca/) Graham cluster of $2.1$GHz Intel $E5-2683$ v4 CPU cores, $32$ cores per node, and each node connected via a $100$ Gb/s network. The results were calculated using the wall clock time taken to simulate $10$ time steps. We ran two tests, using a fixed number of $1$ and $10$ domain decomposition iterations, and we present the scaling results in Fig. 5. The parallel efficiency is very slightly dependent on the number of domain decomposition iterations, and remains above 90% for the number of cores ranging between 32 and 960 (the maximum allocatable without a special permission on the particular cluster). In our opinion, this is an excellent scaling result for an implicit scheme for the incompressible Boussinesq equations.

Acknowledgments

The authors would like to acknowledge the support, under a Discovery Grant, of the National Science and Engineering Research Council of Canada (NSERC).

This research was enabled in part by support provided by Compute Canada (www.computecanada.ca).