An efficient extrapolation multigrid method based on a HOC scheme on nonuniform rectilinear grids for solving 3D anisotropic convection–diffusion problems

Abstract

We develop an efficient multigrid method combined with a high-order compact (HOC) finite difference scheme on nonuniform rectilinear grids for solving 3D diagonal anisotropic convection–diffusion problems with boundary/interior layers. Firstly, we derive a fourth-order compact finite difference scheme to discretize the 3D anisotropic convection–diffusion equation on a rectilinear grid. Then, the resulting large-scale asymmetric linear system of equations is solved by a generalized extrapolation cascadic multigrid (gEXCMG) method based on two novel multigrid (MG) prolongation operators. The highlight of this paper is the application of the quintic Lagrange interpolation and the completed Richardson extrapolation in the design of the new MG prolongation operator on nonuniform rectilinear grids, which can produce a good initial guess (sixth-order approximation to the finite difference solution) for the SSOR-preconditioned biconjugate gradient stabilized (BiCGStab) smoother. In the end, numerical experiments show that the gEXCMG method combined with the HOC scheme can achieve fourth-order accuracy for 3D anisotropic convection–diffusion problems with few smoothing steps on the finest grid. Moreover, the proposed gEXCMG method can offer substantially better efficiency than the state-of-the-art algebraic MG solver, namely, aggregation-based algebraic multigrid (AGMG) method, for large linear systems arising from the discretization of second order elliptic PDEs.

Introduction

The anisotropic convection–diffusion equations have many applications in physical problems such as image processing [1], flows in porous media [2], heat conduction in fusion plasmas [3], transformation thermodynamics [4], and so on. They can also be used to describe some physical phenomena, such as groundwater flow, nanotechnology and transport in biological media [3], [5]. In this study, we consider the steady state solution of some fluid flow and heat transfer problems, which can be described by a three-dimensional (3D) diagonal anisotropic convection–diffusion equation with Dirichlet boundary conditions in the 3D rectangular box $\Omega$ as follows:

$$-a{u}_{xx}-b{u}_{yy}-c{u}_{zz}+p\left(x,y,z\right)u_x+q\left(x,y,z\right)u_y+r\left(x,y,z\right)u_z=f\left(x,y,z\right),\begin{array}{cc}\hfill \hfill & \hfill \left(x,y,z\right)\in \Omega \hfill \end{array}$$$$u\left(x,y,z\right)=g\left(x,y,z\right),\begin{array}{cc}\hfill \hfill & \hfill \left(x,y,z\right)\in \partial \Omega \hfill \end{array}$$

where $a$, $b$ and $c$ denote the diagonal anisotropic coefficients or diffusion coefficients. The coefficients $p\left(x,y,z\right)$, $q\left(x,y,z\right)$, $r\left(x,y,z\right)$ and the forcing function $f\left(x,y,z\right)$, as well as the unknown function $u\left(x,y,z\right)$, are assumed to be continuously differentiable in a rectangular domain $\Omega$ with prescribed boundary values $g\left(x,y,z\right)$ on its boundary $\partial \Omega$. When $a=b=c=1.0$, Eq. (1) reduced to a 3D convection–diffusion equation shown in [6]. When $a=b=c=1.0$, and $p=q=r=0$, Eq. (1) becomes 3D Poisson equation. At present, the traditional numerical solutions used to solve Eq. (1) mainly include finite-difference method, finite-element method, finite-volume method and spectral method [4]. As we all know that solving such a boundary value problem is computationally intensive and challenging especially in three dimensions, even with powerful modern computers, see for example, [4], [7], [8], [9]. With the development of computer technology and the improvement of numerical solution capabilities, people have turned their attention to the HOC difference schemes, which can use fewer grid pedestal points to achieve high-accuracy results [6], [10], [11]. Compared with the traditional high-accuracy numerical solutions, it has advantages of less computational complexity, less sensitivity to elements, and easier handling of boundary conditions [12], [13], [14].

In the past 30 years, the HOC finite difference algorithms have been widely developed to deal with various isotropic and anisotropic convection–diffusion problems [5], [6], [8], [10], [15], [16], [17], [18], [19], [20]. Most existing HOC finite difference schemes are only suitable for uniform grids [13], [21], [22], [23], [24]. If there are not enough grid points in the region of the steep solution gradient, these schemes cannot reach theoretical convergence order for solving the convection–diffusion equations with boundary or internal layers [23], [25], [26], which can be prohibitively computationally expensive. To overcome these issues, efficient HOC finite difference schemes on nonuniform grids are developed [25], [27], [28], [29]. Among these works, one class of nonuniform grid discretization methods is based on the coordinate transformation technique, which transforms the nonuniform grid into a uniform grid [13], [27], [29], [30]. Therefore, we could use the HOC finite difference schemes suitable for the uniform grid to solve these problems [29]. However, these approaches add discrete terms, complicate the derivation, and require the chosen transformation to be explicitly invertible [31], [32]. Another alternative method is directly deriving the HOC finite difference schemes involving no coordinate transformation. For example, Kalita et al. [25] developed a HOC finite difference algorithm to solve the steady 2D convection–diffusion problems with variable coefficients on a rectangular nonuniform grid. Later, it was extended to 3D situations by [6], [33]. These methods have high scale resolutions with a small number of grid points, which can be directly distributed to steep local gradients solution regions or the boundary in the Cartesian coordinate system. However, the HOC differential discretization of 3D convection–diffusion equations results in large and sparse linear equations [34], which constitute the most time-consuming part of the algorithm. Therefore, acceleration algorithms need to be developed to solve these problems.

The MG methods have been presented to efficiently solve the linear algebraic systems generated by HOC finite difference schemes [6], [13], [14], [19], [24], [26], [29], [32], [35], [36], [37]. Since the early HOC finite difference schemes have been developed based on uniform grids or coordinate transformations, most MG methods are implemented on grids with equal mesh sizes [13], [21], [24], [29], [38], [39], [40]. For unequal-meshsize-discretization grids, without coordinate transformation techniques, some specialized multigrid methods based on the partial semi-coarsening strategy have been proposed for the convection–diffusion equation and Poisson’s equation [35], [41]. For general nonuniform rectilinear grids, based on the transformation-free HOC finite difference schemes proposed by [25], Ge and Cao [6], [32] developed an MG V-cycle method to solve the 2D/3D convection–diffusion problems with boundary layers. Shanab et al. [33] adopted the algebraic multigrid (AMG) algorithm to solve HOC finite difference discrete linear systems of 3D convection–diffusion equations on nonuniform rectilinear grids. To our best knowledge, for anisotropic convection–diffusion problems, the standard MG algorithms are not very efficient [37]. Following the idea of [35], Cao et al. [26] proposed a partial semi-coarsening MG approach to solve the convection–diffusion equations with the boundary or interior layers. The computational cost can be greatly reduced without losing accuracy by setting a small number of meshes in non-dominant directions. However, those MG methods have only been tested on problems with a modest number of degrees of freedom [13] and have not been applied to larger-scale problems. Because of the particular high efficiency in solving large-scale linear systems from the elliptic problem [42], [43], [44], [45], [46], the extrapolation cascadic multigrid (EXCMG) method will be adopted to accelerate the solution of the 3D convection–diffusion equations.

The classical MG methods need to cycle between fine and coarse grids. A cascadic multigrid (CMG) method, which is a more straightforward multilevel method without coarse-grid correction, was proposed by Bornemann and Deuflhard [47]. Based on the Richardson extrapolation and CMG method, an EXCMG method was presented to solve the symmetric positive definite (SPD) linear systems arising from the finite element or finite difference discretization of 2D/3D elliptic equations [42], [43], [44]. Recently, a fast cell-centered MG solver based on a finite volume scheme was developed to treat 3D anisotropic diffusion problems with discontinuous coefficients [48]. However, the existing EXCMG method is mainly suitable for solving the SPD linear systems on uniform grids. As far as we know, it has not yet been applied to solve 3D convection–diffusion problems with boundary/interior layers.

This paper proposes a generalized EXCMG (noted as gEXCMG) method combined with HOC finite difference schemes to solve 3D large-scale anisotropic and steady convection–diffusion equations on nonuniform rectilinear grids. Unlike the existing EXCMG methods on uniform grids, the new gEXCMG method can efficiently handle nonuniform grids and problems with large local gradients or boundary/interior layers. For nonuniform grids, based on the sixth-order Lagrange interpolation and completed Richardson extrapolation, we use HOC finite difference solutions at coarse grids to provide a good initial guess for the preconditioned BiCGStab smoother, which can significantly reduce the number of iterations required for each level of grids.

The structure of the paper is as follows. In Section 2, we briefly review the HOC finite difference schemes for 3D diagonal anisotropic convection–diffusion equations on nonuniform rectilinear grids. We explain the principle of the gEXCMG method in detail in Section 3. In Section 4, we verify the proposed gEXCMG algorithm on four numerical experiments. Finally, the conclusions are given in the last section.