A transformation-free HOC scheme and multigrid method for solving the 3D Poisson equation on nonuniform grids

Journal of Computational Physics 234 (2013) 199–216

Contents lists available at SciVerse ScienceDirect

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

A transformation-free HOC scheme and multigrid method for solving the 3D Poisson equation on nonuniform grids Yongbin Ge a,⇑, Fujun Cao b, Jun Zhang c a

Institute of Applied Mathematics and Mechanics, Ningxia University, Yinchuan, Ningxia 750021, China School of Mathematics, Physics and Biological Engineering, Inner Mongolia University of Science and Technology, Baotou, Inner Mongolia 014010, China c Department of Computer Science, University of Kentucky, 773 Anderson Tower, Lexington, KY 40506-0046, USA b

a r t i c l e

i n f o

Article history: Received 28 December 2011 Received in revised form 6 September 2012 Accepted 20 September 2012 Available online 12 October 2012 Keywords: Poisson equations HOC scheme Nonuniform grids Multigrid method Volume law

a b s t r a c t A high-order compact (HOC) difference scheme is proposed to solve the three-dimensional (3D) Poisson equation on nonuniform orthogonal Cartesian grids involving no coordinate transformation from the physical space to the computational space. Theoretically, the proposed scheme has third to fourth-order accuracy; its fourth-order accuracy is achieved under uniform grid settings. Then, a multigrid method is developed to solve the linear system arising from this HOC difference scheme and the corresponding multigrid restriction and interpolation operators are constructed using the volume law. Numerical experiments are conducted to show the computed accuracy of the HOC scheme and the computational efficiency of the multigrid method. Ó 2012 Elsevier Inc. All rights reserved.

1. Introduction In this study, we consider the solution to the three-dimensional (3D) Poisson equation

@2U @2U @2U þ 2 þ 2 ¼ f ðx; y; zÞ @x2 @y @z

ð1Þ

where U is unknown function of the three space variables of x; y and z, f ðx; y; zÞ is the source function, and X a regular domain with @ X is its boundary. The solution Uðx; y; zÞ and f ðx; y; zÞ are assumed to be sufficiently smooth and have the required continuous partial derivatives. The Poisson equation is a basic equation for a lot of fields such as, for instance, in fluid dynamics [1,2], acoustics [3], heat transfer [4], and electromagnetics [5]. With the increasing use of 3D numerical models there is a need for fast algorithms for solving the 3D Poisson equations expressed in finite difference form. Several different methods have been developed, including fast Fourier transform-based (FFT) methods [6,7], cyclic reduction [8,9], the transposed method [10], domain decomposition [11,12], as well as multigrid techniques [13,14]. Among them, multigrid methods combined with high-order compact (HOC) difference schemes are demonstrated to be one of the fastest and most efficient algorithms for solving linear systems arising from the discretization of Poisson equations [15–18] or other elliptic differential equations [19–22]. In [15], Gupta et al. used the nine-point fourth-order compact difference scheme and multigrid V-cycle algorithm to solve the 2D Poisson equation to show the efficiency and performance of this procedure and compared it with the five-point second-order central difference scheme. Zhang [16] employed a fourth-order compact finite difference scheme with multigrid algorithm to solve ⇑ Corresponding author. Tel.: +86 951 2061724; fax: +86 951 2087452. E-mail addresses: [email protected] (Y. Ge), [email protected] (F. Cao), [email protected] (J. Zhang). 0021-9991/$ - see front matter Ó 2012 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.jcp.2012.09.034

200

Y. Ge et al. / Journal of Computational Physics 234 (2013) 199–216

the 3D Poisson equation and compared the influence of different ordering of grid space and different grid transfer operators on the convergence and efficiency of high-order algorithm. Afterwards, Zhang [17] proposed a fourth-order compact discretization scheme with unequal meshsizes and a specialized multigrid method using a partial semi-coarsening strategy and line Gauss–Seidel relaxation to solve the 2D Poisson equation. The two special multigrid methods are found to be very efficient for solving the 2D anisotropic problems. Ge [18] generalized Zhang’s two multigrid methods [17] for solving the 2D Poisson equation to the 3D. We notice that most existing multigrid methods combined with HOC differencing are developed on uniform (usually square or cube) grids [15–22]. However, in many heat and fluid flow problems such as sharp variations of temperature or concentration over a short space, the grids have to be refined over all computational domain if conventional uniform grids are employed. The computation would be prohibitively expensive and wasteful. Under such circumstances, local mesh refinement procedure [13,23–25] or nonuniform grid discretization [26–34] is necessary to make grid points concentrate in regions of sharp variations, thus dramatically reducing the computer storage and computational time. Multigrid local mesh refinement technique was introduced by Brandt [13] in 1977 and then was used to treat Poisson equation with standard discretization scheme [23]. Similar multigrid local mesh refinement procedures were used for solving the CFD computations [24] and the 2D convection diffusion problems with singularity in the source term [25]. Recently, Ge and Cao [34] developed a multigrid method based on the HOC difference scheme on nonuniform grids to solve the 2D convection diffusion problems. The main highlight of this study is that a projection operator and an interpolation operator which is suitable for the HOC scheme on nonuniform grids are presented by the area law, consequently, this is maybe the first time that multigrid technique is seriously considered with the transformation-free HOC scheme on nonuniform grids. In this paper, we are aiming at developing high-order compact difference discretization and multigrid method directly on nonuniform grids for solving the 3D Poisson equation. To the best of our knowledge, there are no reports so far about multigrid method on nonuniform grids with transformation-free HOC difference schemes for the 3D cases. The organization of this paper is as follows. Section 2 presents an HOC finite difference discretization scheme on nonuniform grids for solving the 3D Poisson equation (1); In Section 3, we give a brief introduction to the general philosophy of the multigrid method and propose brand-new restriction and prolongation operators using volume law; In Section 4, numerical experiments are conducted to demonstrate the effectiveness and accuracy of the present method; finally we give concluding remarks in Section 5. 2. HOC scheme on nonuniform grids We consider a regular domain ðx; y; zÞ 2 ½a1 ; a2   ½b1 ; b2   ½c1 ; c2 . The discretization is carried out on a nonuniform 3D grid. Thus, we divide the intervals ½a1 ; a2 , ½b1 ; b2  and ½c1 ; c2  into sub-intervals by the points

a1 ¼ x0 ; x1 ; x2 ; . . . ; xNx 1 ; xNx ¼ a2 ; b1 ¼ y0 ; y1 ; y2 ; . . . ; yNy 1 ; yNy ¼ b2 and

c1 ¼ z0 ; z1 ; z2 ; . . . ; zNz 1 ; zNz ¼ c2 : In the x-direction, letting

hx ¼ ða2  a1 Þ=Nx and the backward and forward step lengths are respectively given by

hbx ¼ xi  xi1 ¼ hbx hx ; hfx ¼ xiþ1  xi ¼ hfx hx ;

1 6 i 6 Nx  1:

And in the y- and z-directions, h^ , hb^ and hf ^ can be defined similarly. ^ denotes x, y or z. To simplify the scheme, furthermore, we define a^ , b^ and c^ as:

a^ ¼ hf ^ hb^ ; b^ ¼ hf ^ þ hb^ ; c^ ¼ hf ^  hb^ ; if and only if

hf ^ ¼ hb^ ¼ 1ðhfx ¼ hbx ; hfy ¼ hby and hfz ¼ hbz Þ the grids turn to be uniform. We use a local coordinate system where the cubic grids are labeled as in Fig. 1. The approximate value of a function Uðx; y; zÞ at an interior grid point ðxi ; yj ; zk Þ is denoted by U0 and the approximate values of its 26 nearest neighboring grid points are defined by Ul l ¼ 1; 2; . . . ; 26, indexed as shown in Fig. 1. For a sufficiently smooth function Uðx; y; zÞ in the given domain, the Taylor series expansion at point 1 gives 2

U1 ¼ U0 þ hfx hx

3

4

5

2 3 4 5 @ U0 hfx hx @ 2 U0 hfx hx @ 3 U0 hfx hx @ 4 U0 hfx hx @ 5 U0 6 þ þ þ þ Oðh6fx hx Þ: þ @x 2! @x2 3! @x3 4! @x4 5! @x5

ð2Þ

201

Y. Ge et al. / Journal of Computational Physics 234 (2013) 199–216

12

20

5

13

11

14

21

19

22

hfz 2

8

(i, j, k)

0

3

7

1

hbz 9

4

10

24

23

16

hfy

17 15

6

25

hby

18

26

hbx

hfx

Fig. 1. The stencil of 3D nonuniform grids.

Similarly at point 3 2

U3 ¼ U0  hbx hx

3

4

5

@ U0 h2bx hx @ 2 U0 h3bx hx @ 3 U0 h4bx hx @ 4 U0 h5bx hx @ 5 U0 6  þ  þ Oðh6bx hx Þ: þ @x 2! @x2 3! @x3 4! @x4 5! @x5

ð3Þ

Using (2) hbx þ (3) hfx , we have 2

3

@ 2 U0 2 hx @ 3 U0 hx 2 @ 4 U0 hx 2 @ 5 U0 4 ¼ ðhbx U1  bx U0 þ hfx U3 Þ  cx   þ Oðgx hx Þ ðbx  3ax Þ ðbx  2ax Þcx 2 2 3 4 @x 3 @x 12 @x 60 @x5 ax bx hx

ð4Þ

in which, gx ¼ b4x  5ax ðb2x  ax Þ. We define the second-order central difference operator in the x-direction as

d2x U0 ¼

2

ax bx h2x

ðhbx U1  bx U0 þ hfx U3 Þ:

ð5Þ

If hbx ¼ hfx ¼ 1, (5) reduces to the central difference operators on uniform grids. Thus, the second-order derivative in the xdirection can be expressed as 2

3

@ 2 U0 hx @ 3 U0 hx 2 @ 4 U0 hx 2 @ 5 U0 4 2 ¼ d U  c   3 a Þ   2 a Þ c þ Oðgx hx Þ: ðb ðb 0 x x x x x x x @x2 3 @x3 12 @x4 60 @x5

ð6Þ

Similarly, we can give the expressions of the second-order derivatives in the y and z-directions. Consequently, (1) can be discretized by

ðd2x þ d2y þ d2z ÞU0 ¼ f0 þ s0 ; where

ð7Þ

s0 is the truncation error given by

s0 ¼ H1

@ 3 U0 @ 3 U0 @ 3 U0 @ 4 U0 @ 4 U0 @ 4 U0 @ 5 U0 @ 5 U0 @ 5 U0 4 þ K1 þ L1 þ H2 þ K2 þ L2 þ H3 þ K3 þ L3 þ Oðgx hx Þ 3 4 3 3 4 4 5 5 @x @y @z @x @y @z @x @y @z5 4

4

þ Oðgy hy Þ þ Oðgz hz Þ; H1 ¼

1 hx c x ; 3

K1 ¼

1 hy cy ; 3

ð8Þ L1 ¼

1 hz c z ; 3

H2 ¼

1 2 2 h ðb  3ax Þ; 12 x x

K2 ¼

1 2 2 h ðb  3ay Þ; 12 y y

202

Y. Ge et al. / Journal of Computational Physics 234 (2013) 199–216

L2 ¼

1 2 2 h ðb  3az Þ; 12 z z

L3 ¼

1 3 2 h ðb  2az Þcz ; 60 z z

If dropping off

ðd2x

þ

d2y

þ

H3 ¼

1 3 2 h ðb  2ax Þcx ; 60 x x

K3 ¼

1 3 2 h ðb  2ay Þcy ; 60 y y

gy ¼ b4y  5ay ðb2y  ay Þ; gz ¼ b4z  5az ðb2z  az Þ:

s0 in (7), we can get the central difference scheme (CDS) on nonuniform grids

d2z Þ

U0 ¼ f0 :

ð9Þ

According to the definition of d2x ; d2y and d2z , the CDS scheme can be written explicitly as

2

1

þ 2

a x hx

! 2hby 2hfx 2hfy 2hfz 1 2hbx 2hbz þ U0 þ U1 þ U2 þ U3 þ U4 þ U5 þ U6 ¼ f0 : ay h2y az h2z ax bx h2x ay by h2y ax bx h2x ay by h2y az bz h2z az bz h2z 1

ð10Þ

Only 7 grid points are involved. From the definition of s0 , we can see than when cx ¼ cy ¼ cz ¼ 0; i.e., hfx ¼ hbx , hfy ¼ hby and hfz ¼ hbz , (10) is of the second-order accuracy. Otherwise, it is just of the first-order accuracy. To further improve the accuracy, we consider

!   2 @ 3 U0 @ 3 U0 @ 3 U0 @ @ @ @ U0 @ 2 U0 @ 2 U0 þ L1 H1 þ K1 þ L1 ¼ H1 þ K 1 þ þ @x @y @z @x3 @y3 @z3 @x2 @y2 @z2 ! @ 3 U0 @ 3 U0 @ 3 U0 @ 3 U0 @ 3 U0 @ 3 U0  H1 þ H1 þ K1 2 þ K1 þ L1 2 þ L1 2 @x@y2 @x@z2 @x @y @y@z2 @x @z @y @z " #   @f0 @f0 @f0 2 @ U2 @ U0 @ U4 2 þ Oðh ¼ H1 h Þ þ K1 þ L1  H1  b þ h by fy y y @x @y @z @x @x @x ay by h2y " #   2 @ U5 @ U0 @ U6 2 þ Oðhz Þ h  H1  b þ h bz fz z @x @x @x az bz h2z " #   2 @ U1 @ U0 @ U3 2  K1 h Þ  b þ h þ Oðh bx fx x x @y @y @y ax bx h2x " #   2 @ U5 @ U0 @ U6 2 h Þ  K1  b þ h þ Oðh bz fz z z @y @y @y az bz h2z " #   2 @ U1 @ U0 @ U3 2 þ Oðhx Þ hbx  L1  bx þ hfx @z @z @z ax bx h2x " #   2 @ U2 @ U0 @ U4 2 h Þ  b þ h þ Oðh  L1 by fy y y @z @z @z ay by h2y

ð11Þ

and

H2

@ 4 U0 @ 4 U0 @ 4 U0 þ K2 þ L2 ¼ 4 4 @x @y @z4

H2

@2 @2 @2 þ K 2 2 þ L2 2 2 @x @y @z

"  ðH2 þ K 2 Þ

!

@ 2 U0 @ 2 U0 @ 2 U0 þ þ @x2 @y2 @z2

!

@ 4 U0 @ 4 U0 @ 4 U0 þ ðH2 þ L2 Þ 2 2 þ ðK 2 þ L2 Þ 2 2 2 2 @x @y @x @z @y @z

#

@ 2 f0 @ 2 f0 @ 2 f0 þ K 2 2 þ L2 2 2 @x @y @z ! " # 2 @ 2 U2 @ 2 U0 @ 2 U4 2 þ Oðhy Þ  ðH2 þ K 2 Þ hby  by þ hfy @x2 @x2 @x2 ay by h2y ! " # 2 @ 2 U5 @ 2 U0 @ 2 U6 2 þ Oðhz Þ hbz  bz þ hfz  ðH2 þ L2 Þ @x2 @x2 @x2 az bz h2z ! " # 2 @ 2 U5 @ 2 U0 @ 2 U6 2 þ Oðhz Þ hbz  bz þ hfz  ðK 2 þ L2 Þ @y2 @y2 @y2 az bz h2z

¼ H2

ð12Þ

The first and second-order derivatives in (11) and (12) can be approximated by the central differences. Then, combining (7) and (8) with (11) and (12), the 19 point high-order compact (HOC) scheme using nonuniform grids for the 3D Poisson equation (1) can be derived as follows

Y. Ge et al. / Journal of Computational Physics 234 (2013) 199–216 18 X @f0 @f0 @f0 @ 2 f0 @ 2 f0 @ 2 f0 Al Ul ¼ f0 þ H1 þ K1 þ L1 þ H 2 2 þ K 2 2 þ L2 2 @x @y @z @x @y @z l¼0

ð13Þ

The coefficients Al ðl ¼ 0; 1; 2; . . . ; 18Þ are given as:

A0 ¼ 2

1 1 þ þ ax h2x ay h2y az h2z hbx

A1 ¼ 2

ax bx h2x hby

A2 ¼ 2

A3 ¼ 2

H1 by

hbz

þ

þ

 2

A10 ¼

A11 ¼

A12 ¼

K 1 bx

ax bx by h2x hy L1 b x

ax bx bz h2x hz

hfz

L1 b x

az bz h2z

þ

ax bx bz h2x hz

2H1 hby

2H1 hfy

ay bx by hx h2y 2H1 hfy 2 y bx by hx hy

a

2H1 hbz 2 z bz bx hx hz

a

2K 1 hbz

az bz by hy h2z

A14 ¼ 





þ

þ



þ



þ

þ

2K 1 hbz

az bz by hy h2z 

az bz by hy h2z

2 x bx by hx hy

a

2L1 hbx 2 x bx bz hx hz

a

2L1 hby

ay bz by h2y hz þ

þ

ax ay bx by h2x h2y ðK 2 þ L2 Þhbz by

4

ay az by bz h2y h2z

!

ðK 2 þ L2 Þhfz by

4

ay az by bz h2y h2z ;

ax ay bx by h2x h2y

4ðH2 þ K 2 Þhfx hby

þ

ax ay bx by h2x h2y 4ðH2 þ K 2 Þhfx hfy

þ

þ

þ

þ

2L1 hfx

2L1 hfy

ay by bz h2y hz 2L1 hbx

4

4ðH2 þ K 2 Þhbx hby

ax bx bz h2x hz

ax bx bz h2x hz

ðH2 þ K 2 Þhfy bx

þ

ax ay bx by h2x h2y

4ðH2 þ K 2 Þhbx hfy

ax ay bx by h2x h2y 4ðH2 þ L2 Þhbx hbz

4ðK 2 þ L2 Þhby hbz

;

;

ay az by bz h2y h2z

þ

4ðH2 þ L2 Þhfx hbz

ax az bx bz h2x h2z 4ðK 2 þ L2 Þhfy hbz

ay az by bz h2y h2z

4ðH2 þ L2 Þhbx hfz

ax az bx bz h2x h2z

;

;

ax az bx bz h2x h2z

þ

;

;

;

;

ðH2 þ L2 Þhbx bz

þ

þ

þ

þ

! ;

ax az bx bz h2x h2z

þ

þ

ax ay bx by h2x h2y

!

ay bz by h2y hz

2K 1 hfx

ðH2 þ K 2 Þhfx by

!

K 1 bz

ax bx by h2x hy

ax ay bx by h2x h2y

4

az bz bx hx h2z

ay bz by h2y hz

ðH2 þ K 2 Þhby bx

4

!

H1 bz

L1 by

ax ay bx by h2x h2y

!

az bz by hy h2z

L1 by

ðH2 þ K 2 Þhbx by

4

K 1 bz

þ

2K 1 hbx

!

az bz bx hx h2z

ax bx by h2x hy



! K 2 þ L2 H2 þ L2 þ þ ; ax ay h2x h2y ay az h2y h2z ax az h2x h2z H2 þ K 2

H1 bz

2K 1 hfx

þ

2H1 hbz

2H1 hfz

þ4

ax bx by h2x hy

az bz bx hx h2z

az bz bx hx h2z

!

2K 1 hlx

þ

ay bx by hx h2y

A13 ¼ 

A15 ¼

ay bx by hx h2y

az bz hz

2H1 hby

A9 ¼ 

K 1 bx

hfx

ay bx by hx h2y

A8 ¼ 

ay bx by hx h2y

ax bx by h2x hy

ay by h2y

A6 ¼ 2

H1 by

ay by hy

hfy

A5 ¼ 2



 2

ax bx h2x

A4 ¼ 2

A7 ¼

1

203

ðK 2 þ L2 Þhby bz

! ;

ay az by bz h2y h2z

ðH2 þ L2 Þhfx bz

! ;

ax az bx bz h2x h2z ðK 2 þ L2 Þhfy bz

! ;

ay az by bz h2y h2z

ðH2 þ L2 Þhbz bx

! ;

ax az bx bz h2x h2z ðH2 þ L2 Þhfz bx

ax az bx bz h2x h2z

! ;

204

Y. Ge et al. / Journal of Computational Physics 234 (2013) 199–216

A16 ¼

2K 1 hfz 2 z bz by hy hz

a

A17 ¼ 

A18 ¼ 



2H1 hfz 2 z bz bx hx hz

a

2K 1 hfz

az bz by hy h2z

2L1 hby 2 y by bz hy hz

a





þ

2L1 hfx 2 x bx bz hx hz

a

2L1 hfy

ay by bz h2y hz

4ðK 2 þ L2 Þhby hfz

ay az by bz h2y h2z þ

þ

;

4ðH2 þ L2 Þhfx hfz

ax az bx bz h2x h2z 4ðK 2 þ L2 Þhfy hfz

ay az by bz h2y h2z

;

:

From expressions of s0 , it is easy to know that the accuracy of the scheme is of the third- to fourth-order; its fourth-order accuracy is achieved on uniform grids when c^ ¼ 0; i.e., hb^ ¼ hf ^ ¼ 1. Under this condition, the HOC scheme reduces to the fourth-order compact scheme proposed by Wang et al. [35], which is also derived by Ge [18]. Otherwise, if hb^ – 1 or hf ^ – 1, at least third-order accuracy is achieved on nonuniform grids. 3. Multigrid method Multigrid method is one of the fastest and most efficient iterative methods for solving a wide class of partial differential equations. This method offers convergence rate independent of the grid size and is very effective for solving large scale sparse linear systems which are obtained by discretizing elliptic problems [13–22,25,34,36,37]. The essential principle of multigrid method is to approximate the smooth (long wavelength) part of the error on coarser grids. The non-smooth or rough part is reduced with a small number (independent of meshsize h) of iterations with a basic iterative method (e.g., Gauss–Seidel or Jacobi methods, etc.) on the fine grid. One iteration of a simple multigrid cycling algorithm consists of smoothing the error using a relaxation technique, restricting the residuals to the coarse grid, solving an approximation to the smooth error equation on the coarse grid, interpolating the coarse grid error correction back to the fine grid and finally adding the error correction into the approximation. Hence, restriction, interpolation and relaxation operators are three important components in multigrid method and only when they are suitably chosen and match well with each other, can the efficiency of multigrid method be optimal. For more details about multigrid method, readers are referred to [13,14,36,37] and the references therein. For solving the 2D and 3D Poisson equations with fourth-order compact discretization schemes on uniform grids, some specific multigrid methods are employed [15–18]. A standard bilinear interpolation operator and a full weighting restriction operator are used as the inter grid transfer operators. However, on nonuniform grids those restriction and interpolation operators developed on uniform grids [36,37] are not suitable any longer. If we still use them, the inexactness of them will cause a significant decrease in the efficiency of the standard multigrid methods and the convergence rate will deteriorate considerably. Therefore, we will develop a new restriction operator and interpolation operator which are suitable for the HOC scheme on nonuniform grids. 3.1. Restriction operator The principle of constructing a restriction operator is to evaluate the residuals on coarse grids with the help of the residuals on fine grids. Liu [38] developed a so called area law for the restriction of residual in multigrid method in 2D. Although finite volume schemes on uniform grids are used therein, the methodology can be generalized to construct both restriction and interpolation operators for finite difference schemes on nonuniform grids. For 3D cases, we call the methodology volume law. For each grid point on the coarse grid levels, besides itself, there are correspondingly 26 fine grid points surrounding it. Since the reference grid point and the surrounding grid points on the fine grids have contributions of different degree to it on the coarse grids, a kind of full-weighting restriction operator based on the volume law for nonuniform grids is constructed. For convenience, we mark these 27 grid points with order as in Fig. 1. The volumes concerning with each grid points on nonuniform grids are bounded by half meshsize lines around the reference grid points. For instance, the shadow in Fig. 2(a) gives the definition of V 0 and the shadow in Fig. 2(b) gives the definition of V 1 and V 3 . The volumes of the other grid points can be defined similarly. V is the total volume of V i (i ¼ 0; . . . ; 26). The key point of deriving the full weighting restriction operator is to evaluate the weighting coefficients of the residuals of each grid points. The reference point ði; j; kÞ on the fine grids has the most contributions to it on the coarse grids, so the corresponding weighting coefficient is evaluated by V 0 =V. At the same time, we notice that the grid points near the reference point ði; j; kÞ (with smaller volumes) have much more contributions than those far away from it (with bigger volumes); i.e., we can use the volumes in opposite direction proportional to the total volume V as the weighting coefficient of the present grid point. For instance, the weighting coefficient of the point ði þ 1; j; kÞ is given by V 3 =V and the point ði  1; j; kÞ by V 1 =V, etc. Hence, suppose that ri;j;k is the residual at the fine grid point ði; j; kÞ and  It is easy to know that i ¼ 2i, j ¼ 2j, k ¼ 2k. The full weighting ri;j;k is the corresponding residual at the coarse grid point ði; j; kÞ. restriction operator on nonuniform grids can be explicitly written as

Y. Ge et al. / Journal of Computational Physics 234 (2013) 199–216

205

(b)

(a)

Fig. 2. The restriction operator on nonuniform grids.

ri;j;k ¼

1 ðV 0 r i;j;k þ V 1 ri1;j;k þ V 2 ri;j1;k þ V 3 r iþ1;j;k þ V 4 r i;jþ1;k þ V 5 r i;j;k1 þ V 6 ri;j;kþ1 þ V 7 ri1;j1;k þ V 8 riþ1;j1;k þ V 9 r iþ1;jþ1;k V þ V 10 r i1;jþ1;k þ V 11 r i1;j;k1 þ V 12 r i;j1;k1 þ V 13 r iþ1;j;k1 þ V 14 r i;jþ1;k1 þ V 15 ri1;j;kþ1 þ V 16 ri;j1;kþ1 þ V 17 riþ1;j;kþ1 þ V 18 r i;jþ1;kþ1 þ V 19 r i1;j1;k1 þ V 20 r iþ1;j1;k1 þ V 21 riþ1;jþ1;k1 þ V 22 r i1;jþ1;k1 þ V 23 r i1;j1;kþ1 þ V 24 r iþ1;j1;kþ1 þ V 25 r iþ1;jþ1;kþ1 þ V 26 ri1;jþ1;kþ1 Þ:

in which,

V ¼ ðhfx þ hbx Þ  ðhfy þ hby Þ  ðhfz þ hbz Þ; V0 ¼

1 ðhfx þ hbx Þ  ðhfy þ hby Þ  ðhfz þ hbz Þ; 8

V2 ¼

1 hfy  ðhfx þ hbx Þ  ðhfz þ hbz Þ; 8

V3 ¼

1 hbx  ðhfy þ hby Þ  ðhfz þ hbz Þ; 8

V4 ¼

1 hby  ðhfx þ hbx Þ  ðhfz þ hbz Þ; 8

V5 ¼

1 hfz  ðhfx þ hbx Þ  ðhfy þ hby Þ; 8

V6 ¼

1 hbz  ðhfx þ hbx Þ  ðhfy þ hby Þ; 8

V7 ¼

1 hfx  hfy  ðhfz þ hbz Þ; 8

V8 ¼

1 hbx  hfy  ðhfz þ hbz Þ; 8

V9 ¼

V1 ¼

1 hfx  ðhfy þ hby Þ  ðhfz þ hbz Þ; 8

1 hbx  hby  ðhfz þ hbz Þ; 8

V 10 ¼

1 hfx  hby  ðhfz þ hbz Þ; 8

V 11 ¼

1 hfx  hfz  ðhfy þ hby Þ; 8

V 12 ¼

1 hfy  hfz  ðhfx þ hbx Þ; 8

V 13 ¼

1 hbx  hfz  ðhfy þ hby Þ; 8

V 14 ¼

1 hby  hfz  ðhfx þ hbx Þ; 8

V 15 ¼

1 hfx  hbz  ðhfy þ hby Þ; 8

V 16 ¼

1 hfy  hbz  ðhfx þ hbx Þ; 8

V 17 ¼

1 hbx  hbz  ðhfy þ hby Þ; 8

206

Y. Ge et al. / Journal of Computational Physics 234 (2013) 199–216

V 18 ¼

1 hby  hbz  ðhfx þ hbx Þ; 8

V 20 ¼

1 hbx  hfy  hfz ; 8

V 21 ¼

1 hbx  hby  hfz ; 8

V 22 ¼

1 hfx  hby  hfz ; 8

V 23 ¼

1 hfx  hfy  hbz ; 8

V 24 ¼

1 hbx  hfy  hbz ; 8

V 25 ¼

1 hbx  hby  hbz ; 8

V 26 ¼

1 hfx  hby  hbz : 8

V 19 ¼

1 hfx  hfy  hfz ; 8

If the grid spacing reduces to be equal meshsize, then the total volume V is equally divided into sixty-four small parts by the grid lines and the half meshsize lines. Suppose the volume of each small part is V, we can get that V 0 ¼ 8V, V 1 ¼ V 2 ¼    ¼ V 6 ¼ 4V, V 7 ¼ V 8 ¼    ¼ V 18 ¼ 2V and V 19 ¼ V 20 ¼    ¼ V 26 ¼ V. Under this condition, the restriction operator reduces to the conventional full weighting operators on equal-meshsize-discretization grids [37]

ri;j;k ¼

1 ½8r i;j;k þ 4ðri1;j;k þ riþ1;j;k þ r i;j1;k þ r i;jþ1;k þ r i;j;k1 þ ri;j;kþ1 Þ þ 2ðr iþ1;jþ1;k þ ri1;jþ1;k þ r iþ1;j1;k þ r i1;j1;k 64 þ r i;jþ1;k1 þ ri;jþ1;kþ1 þ ri;j1;k1 þ r i;j1;kþ1 þ r iþ1;j;kþ1 þ r iþ1;j;k1 þ r i1;j;k1 þ r i1;j;kþ1 Þ þ ðr iþ1;jþ1;k1 þ r i1;jþ1;k1 þ r iþ1;j1;k1 þ r i1;j1;k1 þ r iþ1;jþ1;kþ1 þ r i1;jþ1;kþ1 þ r iþ1;j1;kþ1 þ r i1;j1;kþ1 Þ:

3.2. Interpolation operator For interpolation operator, we use a similar strategy to construct it. We notice that the grid points on coarse grid level are simultaneously the grid points on fine grid level when grids are transferred from the coarse grids to the fine grids. Thus, we directly transfer these grid points from the coarse grid to the fine grid, see Fig. 3(a). The interpolation operator is explicitly written as

ri;j;k ¼ ri;j;k : Then, for those grid points on the fine grid (the small black dots) between two coarse grids (the big black dots) we can interpolate them with their two neighbor grid points on coarse grids, see Fig. 3(b). The error corrections along x-, y- and zdirections are respectively interpolated by the following formula

ri1;j;k ¼

1 ðhfxri1;j;k þ hbxri;j;k Þ; hfx þ hbx

ri;j1;k ¼

1 ðhfyri;j1;k þ hbyr i;j;k Þ; hfy þ hby

ri;j;k1 ¼

1 ðhfzr i;j;k1 þ hbzri;j;k Þ:  hfz þ hbz

Thirdly, for those central grid points in the even planes (the small black dots), see Fig. 3(c), we use four grid points around them (the big black dots) on the coarse grids to interpolate as following

ri1;j1;k ¼

1 ðS1xyri1;j1;k þ S2xyri;j1;k þ S3xyri;j;k þ S4xyri1;j;k Þ; Sxy

ri;j1;k1 ¼

1 ðS1yzri;j1;k1 þ S2yzri;j;k1 þ S3yzr i;j;k þ S4yxr i;j1;k Þ;   Syz

ri1;j;k1 ¼

1 ðS1xzri1;j;k1 þ S2xzri;j;k1 þ S3xzri;j;k þ S4xyri1;j;k Þ   Sxz

in which,

Sxy ¼ ðhfx þ hbx Þ  ðhfy þ hby Þ; Syz ¼ ðhfy þ hby Þ  ðhfz þ hbz Þ; Sxz ¼ ðhfx þ hbx Þ  ðhfz þ hbz Þ; S1xy ¼ hfx  hfy ; S2xy ¼ hbx  hfy ; S3xy ¼ hbx  hby ; S4xy ¼ hfx  hby ;

Y. Ge et al. / Journal of Computational Physics 234 (2013) 199–216

(a)

(b)

(c)

(d)

207

Fig. 3. The interpolation operator on nonuniform grids.

S1yz ¼ hfy  hfz ; S2yz ¼ hby  hfz ; S3yz ¼ hby  hbz ; S4zy ¼ hfy  hbz ; S1xz ¼ hfx  hfz ;

S2xz ¼ hbx  hfz ;

S3xz ¼ hbx  hbz ;

S4xz ¼ hfx  hbz :

Finally, for the central grid points in the odd planes (the small black dots), see Fig. 3(d), we use the eight grid points around them (the big black dots) on the coarse grids to make average as following

ri1;j1;k1 ¼

1 ~ ~ 2r   þ V~ 3r   þ V ~ 4r   þ V ~ 5r   þ V~ 6r   þ V ~ 7r   þ V~ 8r   Þ; ðV r   þ V i;j1;k1 i;j;k1 i1;j;k1 i1;j1;k i;j1;k i;j;k i1;j;k ~ 1 i1;j1;k1 V

~ i (i ¼ 1; . . . ; 8) is the volume formed by all the grid lines around the reference point ði  1; j  1; k  1Þ, respectively. V ~ where V ~ i (i ¼ 1; . . . ; 8) and they are calculated by following formula is the total volume of V

~ ¼ ðhfx þ hbx Þ  ðhfy þ hby Þ  ðhfz þ hbz Þ; V

208

Y. Ge et al. / Journal of Computational Physics 234 (2013) 199–216

~ 1 ¼ hfx  hfy  hfz ; V

~ 2 ¼ hbx  hfy  hfz ; V

~ 5 ¼ hfx  hfy  hbz ; V

~ 6 ¼ hbx  hfy  hbz ; V

~ 3 ¼ hbx  hby  hfz ; V

~ 4 ¼ hfx  hby  hfz ; V

~ 7 ¼ hbx  hby  hbz ; V

~ 8 ¼ hfx  hby  hbz : V

If the grid spacing reduces to equal meshsize, the interpolation operator reduces to the conventional trilinear interpolation operator on equal-meshsize-discretization grids [37]

ri;j;k ¼ ri;j;k ; ri1;j;k ¼

1 ðr þ ri;j;k Þ; 2 i1;j;k

ri;j1;k ¼

1 ðr  þ r i;j;k Þ; 2 i;j1;k

ri;j;k1 ¼

1 ðr  þ r i;j;k Þ; 2 i;j;k1

ri1;j1;k ¼

1 ðr  þ ri;j1;k þ ri;j;k þ ri1;j;k Þ; 4 i1;j1;k

ri;j1;k1 ¼

1 þ r i;j;k þ r i;j1;k Þ; ðr   þ r i;j;k1  4 i;j1;k1

ri1;j;k1 ¼

1 ðr  þ ri;j;k1 þ ri;j;k þ ri1;j;k Þ;  4 i1;j;k1

ri1;j1;k1 ¼

1 ðr   þ ri;j1;k1 þ ri;j;k1 þ ri1;j;k1 þ ri1;j1;k þ ri;j1;k þ ri;j;k þ ri1;j;k Þ:    8 i1;j1;k1

3.3. Relaxation operator Relaxation operator (the smoother) is another important operator in the multigrid method. Its role is not to remove the errors, but to damp the high frequency components of the errors on the current grid while leaving the low frequency components to be removed by the coarser grids. For simple, isotropic problems, even simple smoothers (e.g., pointwise Gauss– Seidel relaxations [15,16]) efficiently damp a smooth error in all coordinate directions. But for anisotropic, boundary layers or local singularities problems, line Gauss–Seidel relaxation [17,18] and alternating line Gauss–Seidel relaxation [33,34,39,40] are verified to be more robust smoothers. As for 3D cases, the point Gauss–Seidel relaxation in lexicographical ordering [16], red and black ordering [16,22], four color ordering [21,22] and the alternating direction plane relaxation [41,42] are used by several authors. We note that four

z

R

B

R

B

R

B

R

B

G

O

G

O G

O

G

O

R

B

R

B

B

R

B

top plane

x R

y O

G

O

G

O

G

O

G

B

R

B

R

B

R

B

R

O

G

O

G

O

G

O

G

R

B

R

B

R

B

R

B

G

O

G

O G

O

G

O

R

B

R

B

B

R

B

R

Fig. 4. Decoupling the 3D grid-points with four colors.

middle plane

bottom plane

209

Y. Ge et al. / Journal of Computational Physics 234 (2013) 199–216

colors are sufficient to completely decouple the 3D grid points with 19 point compact scheme, so four-color Gauss–Seidel relaxation exhibits good performance for the computations with the HOC scheme. Fig. 4 depicts the decoupling of the 3D grid-points with four colors, in which grid points colored with red (R), black (B), green (G) and orange (O) colors are used. Another efficient smoother is the plane relaxation [41,42], which update all grid points in a plane simultaneously. If the planes are swept in all directions (e.g., in the xy-, xz-, and yz-planes) sequentially, the resulting error is smooth in all directions. The increased cost of this method lies, obviously, in the increased cost of the smoother. In this paper, we use these four relaxations to smooth the residuals on each coarse grid. For the alternating direction plane relaxation, we perform one sweep of the xy-plane Gauss–Seidel relaxation along the z-coordinate direction first, followed by one sweep of the yz-plane GaussSeidel relaxation and xz-plane Gauss–Seidel relaxation along the x-coordinate direction and y-coordinate direction, respectively. On each 2D grid plane, we use one line Gauss–Seidel relaxation [17,18]. Since all the grid points are computed three times for each time alternating direction plane relaxation but only one time for each time point Gauss–Seidel relaxation, red–black Gauss–Seidel relaxation or four-color Gauss–Seidel relaxation, we use multigrid V(1, 1) cycle for the alternating direction plane relaxation (Plane GS) but multigrid V(3, 3) cycle for point Gauss–Seidel relaxation (PGS), red–black Gauss– Seidel relaxation (RBGS) and four-color Gauss–Seidel relaxation (four-color GS) to keep approximately equal computational cost for each time multigrid V-cycle. In the next section, we will study smoothing effect of different relaxation operators by numerical experiments. 4. Numerical experiments Three test problems are chosen to numerically validate the accuracy and effectiveness of the present method. The multigrid V-cycle process is used and the computation is started with zero initial guesses and is terminated when the Euclidean norm (2-norm) of the initial residual vector on the finest grids is reduced by 1010. The code is written in Fortran 77 programming language with double precision arithmetic. All computations are run on a personal computer with an Intel 2.4 GHz CPU and 2 GB memory. 4.1. Test problem 1 Consider the following Poisson problem with a source term f

Uxx þ Uyy þ Uzz ¼ f ðx; y; zÞ; f ðx; y; zÞ ¼ 

0 < x; y; z < 2;

px py pz 3p2 sin sin : sin 2 2 2 4

On all the boundaries, U is supposed to be zero. Its analytic solution is

Uðx; y; zÞ ¼ sin

px 2

sin

py 2

sin

pz 2

:

We notice that the exact solution of this problem does not exhibit any sharp variations, so nonuniform grids are not necessary. We just use uniform grids for this test problem to validate the computed accuracy that is achievable by the HOC scheme and the effectiveness of multigrid method introduced in this paper. For comparison, the results of the CDS scheme (10) are also given. The errors reported are the maximum absolute errors over the entire discretized grid points on the finest grid. The accuracy order of a difference scheme can be evaluated by the following formula

Order ¼ log2

ErrorðN1 Þ ; ErrorðN2 Þ

where ErrorðN 1 Þ and ErrorðN 2 Þ are the maximum absolute errors estimated for two different grids with N 1 þ 1 and N 2 þ 1 points in each directions and N 2 is twice N 1 . Firstly, we use different meshsizes from 83 to 1283 grids to evaluate the computed accuracy order. Table 1 gives the maximum absolute errors and the estimated accuracy order for Test problem 1. Since this is an isotropic problem, a very accurate

Table 1 Maximum absolute errors and the accuracy orders of the two schemes for Test problem 1. Grids

CDS Error

83

1.30(2)

HOC Order

Error

Order

7.47(4)

163

3.22(3)

2.01

4.57(5)

4.03

323

8.04(4)

2.00

2.84(6)

4.01

643

2.01(4)

2.00

1.77(7)

4.00

1283

5.02(5)

2.00

1.11(8)

4.00

210

Y. Ge et al. / Journal of Computational Physics 234 (2013) 199–216

Table 2 The number of multigrid V-cycles and CPU time with two schemes and different relaxation methods for Test problem 1. Grids

CDS

Red–black GS

Four-color GS

Plane GS

CPU

Iter.

CPU

Iter.

CPU

Iter.

CPU

3

11

0.094

10

0.078

10

0.079

10

0.109

323

11

0.391

11

0.390

10

0.406

10

0.453

643

12

2.500

9

1.813

10

2.516

10

4.328

163

8

0.078

8

0.078

7

0.063

8

0.094

323

9

0.438

8

0.485

8

0.406

8

0.578

643

9

2.703

8

3.531

8

2.688

8

5.157

16

HOC

Point GS Iter.

solution is obtained on uniform grids with the HOC scheme. Besides, we can see that both the HOC scheme and the CDS scheme obtain their optimal accuracy order on uniform grids and the HOC scheme gets more accurate solution than the CDS scheme. Table 2 lists the number of multigrid V-cycles and the corresponding CPU time in seconds for solving Test problem 1 on the 163 , 323 and 643 grids. We can see that for this problem the multigrid method is very efficient and all the smoothers work well. The computed results show that for the HOC scheme the four-color Gauss–Seidel relaxation and the Point Gauss–Seidel relaxation are more cost -effective (less CPU time) than the red–black Gauss–Seidel relaxation and the alternating plane Gauss–Seidel relaxation. For the CDS scheme, the red–black Gauss–Seidel relaxation shows more cost-effective than the other three smoothers. For both the HOC scheme and the CDS scheme the alternating plane Gauss–Seidel relaxation cost the most CPU time, so it is the most expensive. So, we can say that pointwise relaxations are better choices than plane relaxations for isotropic problems on uniform grids. 4.2. Test problem 2 Consider the following differential equation in the presence of source term

Uxx þ Uyy þ Uzz ¼ f ðx; y; zÞ;

0 < x; y; z < 1;

f ðx; y; zÞ ¼ 104 e100x ½yð1  yÞ  x  zð1  zÞ þ 200e100x for which the analytic solution is

Uðx; y; zÞ ¼ e100x ½yð1  yÞ  x  zð1  zÞ: The boundary conditions are given by the analytic solution. This problem has a steep boundary layer along x ¼ 0. Therefore, a nonuniform grid along the x-direction with grid clustering near x ¼ 0 is used and uniform grid along the y- and zdirections is kept by the following stretching function (see [32])

xi ¼

  i k pi ; þ sin Nx p Nx

yj ¼

j ; Ny

zk ¼

k ; Nz

where k is stretching parameter controlling the density of grid points in the x direction. k < 0 indicates that more grid points are clustered to the boundary x ¼ 0 and similarly to the boundary x ¼ 1 for k > 0. If k ¼ 0 grids turn to be uniform. When the mesh is 323 and k ¼ 0:8, the grids distribution in the xy-plane is given as in Fig. 5. Table 3 gives the maximum absolute errors and the estimated accuracy orders with different stretching parameter k for Test problem 2. k changes from 0.0 to 0.9. We can see that on uniform grid (k ¼ 0:0), the computed results are not very accurate (big errors) and the convergence orders cannot reach four for the HOC scheme or two for the CDS scheme as theoretically predicted. On the contrary, on the nonuniform grid with the decreasing of stretching parameter k, more and more grid points are clustered into the boundary layers, consequently, more and more accurate solution and convergence orders are obtained for both difference schemes. We notice that when k ¼ 0:8, the most accurate solution (the smallest error) is obtained with the HOC scheme. But when k continuously decreases to 0:9, the accuracy decreases. This is not surprising because excessively putting more grid points into the boundary layer area will inevitably cause lack of the grid points in the other regions in the computational domain. Fig. 6 shows the contours of exact solution in the plane y ¼ 0:5 (a), computed solution by the CDS scheme on uniform grid (b), by the HOC scheme on uniform grid (c), and by the HOC scheme on nonuniform grid with k ¼ 0:8 (d). We can see that both the CDS scheme and the HOC scheme on uniform grid get distorted solutions near boundary layers, but the HOC scheme on nonuniform grid with k ¼ 0:8 captures very well with the exact solution curve and almost indistinguishable from the exact solution. Table 4 lists the number of multigrid V-cycles with different relaxation methods and the corresponding CPU time in seconds for different stretching parameter k for Test problem 2. We find that on uniform grids or on nonuniform grids with the grid stretching parameter k being big; e.g., k ¼ 0:0; 0:2; 0:4 (small stretching ratio for grids), all smoothers remove the er-

211

Y. Ge et al. / Journal of Computational Physics 234 (2013) 199–216

Fig. 5. Nonuniform grids distribution in the xy-plane (323 , k ¼ 0:8).

Table 3 Maximum absolute errors and the accuracy orders of the two schemes for Test problem 2. Scheme

k

323

643

Order

1283

Order

CDS

0.0 0.2 0.4 0.6 0.8 0.9

6.49(2) 5.11(2) 3.26(2) 1.59(2) 3.41(3) 7.71(4)

2.34(2) 1.62(2) 9.42(3) 4.24(3) 8.72(4) 1.92(4)

1.47 1.67 1.79 1.91 1.97 2.01

6.73(3) 4.34(3) 2.47(3) 1.09(3) 2.19(4) 4.80(5)

1.79 1.90 1.93 1.96 1.99 2.00

HOC

0.0 0.2 0.4 0.6 0.8 0.9

2.04(2) 1.07(2) 3.95(3) 8.05(4) 4.96(5) 1.43(4)

2.04(3) 9.08(4) 2.95(4) 5.50(5) 3.08(6) 8.98(6)

3.32 3.56 3.74 3.87 4.00 4.00

1.50(4) 6.18(5) 1.95(5) 3.54(6) 1.92(7) 5.60(7)

3.77 3.88 3.92 3.96 4.00 4.00

rors efficiently. But for small k; e.g., k ¼ 0:6; 0:8; 0:9 (big stretching ratio for grids), more multigrid V-cycle numbers are needed with pointwise smoothers while no significant increase for the alternating plane Gauss–Seidel relaxation for the CDS scheme. The same multigrid V-cycle numbers are taken for the HOC scheme under such circumstance. We notice that the multigrid method with the four-color Gauss–Seidel relaxation is divergent for k ¼ 0:8 and 0:9 for the CDS scheme. One possible explanation is that the four-color Gauss–Seidel relaxation based on the nonuniform CDS scheme cannot efficiently damp the smooth errors on the coarse grids, which leads to the correction of the errors on the fine grids to be inaccurate. The inaccuracy becomes more and more serious with multigrid V-cycles, which causes the multigrid method to diverge at last. For all k, the alternating plane Gauss–Seidel relaxation is the most efficient smoother with the least multigrid V-cycle numbers. As far as the CPU time is concerned, the red–black Gauss–Seidel relaxation and the four-color Gauss–Seidel relaxation take less CPU time and thus is more cost -effective for uniform grids or nonuniform grids with big k. But this advantage is lost for small k. In this case, the alternating plane Gauss–Seidel relaxation takes the least CPU time with the least multigrid iteration numbers. 4.3. Test problem 3 Consider the following differential equation in the presence of source term

uxx þ uyy þ uzz ¼ f ðx; y; zÞ; for which the analytic solution is

0 < x; y; z < 1

212

Y. Ge et al. / Journal of Computational Physics 234 (2013) 199–216

0.25

0.25

0.2

0.2

0.1

solution

Solution

0.15

0.15 0.1

0.05

0.05

0 0

0

0.2

0.2

0.2 0.6

0.6 0.8

0.8

0.2 0.4

0.4

0.4

z

0 0

0 0.4

z

x

0.6

0.6 0.8

0.8

x

1 1

1 1

(b)

(a)

0.25

0.25

0.2

0.2

0.1

solution

solution

0.15

0.15 0.1

0.05

0.05

0 0

0

z

0.2

0.4

0.4 0.6

0.6 0.8

0.8

0 0

0

0.2

0.2

0.2 0.4

x

z

0.4 0.6

0.6 0.8

1 1

0.8 1

x

1

(d)

(c)

Fig. 6. Solutions of Problem 2: (a) exact, (b) CDS scheme on uniform grid, (c) HOC scheme on uniform grid, and (d) HOC scheme on nonuniform grid ðk ¼ 0:8Þ.

Table 4 The number of multigrid V-cycles and CPU time with two schemes and different relaxation methods with 323 for Test problem 2. Scheme

k

Point GS

Red–black GS

Four-color GS

Plane GS

Iter.

CPU

Iter.

CPU

Iter.

CPU

Iter.

CPU

CDS

0.0 0.2 0.4 0.6 0.8 0.9

9 10 11 11 21 75

0.344 0.437 0.484 0.515 0.954 3.350

7 8 10 11 18 22

0.296 0.360 0.453 0.500 0.844 1.031

7 7 10 21 – –

0.296 0.343 0.469 0.969 – –

9 8 7 8 9 10

0.579 0.625 0.500 0.547 0.719 0.797

HOC

0.0 0.2 0.4 0.6 0.8 0.9

9 8 8 10 11 15

0.906 0.899 0.899 1.125 1.224 1.578

8 8 8 10 11 15

0.891 0.894 0.894 1.109 1.219 1.562

8 8 8 10 11 16

0.857 0.859 0.859 1.096 1.218 1.617

7 7 7 7 7 7

0.875 0.897 0.897 0.897 0.897 0.897

‘–’ indicates lack of convergence.

213

Y. Ge et al. / Journal of Computational Physics 234 (2013) 199–216

uðx; y; zÞ ¼

½1  e100ðx1Þ ½1  e100ðy1Þ ½1  e100ðz1Þ  ð1  e100 Þ3

:

The source term is determined by the exact solution, which has boundary layers along x ¼ 1, y ¼ 1 and z ¼ 1. Therefore, a nonuniform grid along all three directions with clustering near x ¼ 1, y ¼ 1 and z ¼ 1 is used by the following stretching functions

xi ¼

i k pi þ sinð Þ; Nx p Nx

yj ¼

  j k pj ; þ sin Ny p Ny

zk ¼

  k k pk : þ sin Nz p Nz

When k is more close to 1, more grid points are clustered near x ¼ 1, y ¼ 1 and z ¼ 1. When the mesh is 323 and k ¼ 0:8, the grids distribution in the xy -plane is given as in Fig. 7. Table 5 gives the maximum absolute errors and the estimated accuracy orders with different stretching parameter k for Test problem 3. k changes from 0.0 to 0.9. We can see that on uniform grid (k ¼ 0:0), the convergence orders cannot reach four for the HOC scheme or two for the CDS scheme. On the contrary, on the nonuniform grids with the increasing of stretching parameter k, more and more grid points are clustered into the boundary layers, more and more accurate solution is obtained from the two difference schemes. And convergence orders increase continuously with the increase of k. We notice that when k ¼ 0:8, the most accurate solution is achieved with the HOC scheme. But when k continuously increases to 0.9, the accuracy decreases. Fig. 8 depicts the contours of exact solution in the plane y ¼ 0:6875 (a), computed solution by the CDS scheme on uniform grid (b), by the HOC scheme on uniform grid (c), and by the HOC scheme on nonuniform grid with k ¼ 0:8 (d). We can see that the CDS scheme gets very poor solution. The solution distribution of the HOC scheme on uniform

Fig. 7. Nonuniform grids distribution in the xy-plane (323 , k ¼ 0:8).

Table 5 Maximum absolute errors and the accuracy orders of the two schemes for Test problem 3. Scheme

k

323

643

Order

1283

Order

CDS

0.0 0.2 0.4 0.6 0.8 0.9

5.28(1) 3.89(1) 2.46(1) 1.18(1) 2.70(2) 4.26(3)

1.80(1) 1.20(1) 6.96(2) 3.13(2) 6.84(3) 1.07(3)

1.55 1.70 1.82 1.91 1.98 1.99

4.92(2) 3.18(2) 1.80(2) 7.93(3) 1.72(3) 2.67(4)

1.87 1.92 1.95 1.98 1.99 2.00

HOC

0.0 0.2 0.4 0.6 0.8 0.9

1.45(1) 7.15(2) 2.61(2) 5.18(3) 5.06(4) 1.26(3)

1.38(2) 5.92(3) 1.92(3) 3.51(4) 3.12(5) 7.87(5)

3.39 3.59 3.76 3.88 4.02 4.00

9.66(4) 3.99(4) 1.25(4) 2.24(5) 1.94(6) 4.91(6)

3.84 3.89 3.94 3.97 4.01 4.00

214

Y. Ge et al. / Journal of Computational Physics 234 (2013) 199–216

1

1

0.8

0.4

0.8

Solution

Solution

0.6

0.6 0.4

0.2

0.2

0 0

0

z

0 0

0

0.5

0.5

z

x

0.5

0.5

1 1

11

(a)

(b)

x

1

1

0.8

0.8

0.4

Solution

Solution

0.6

0.6

0.4

0.2 0 0

0

z

0.5

0.5

0.2

0 0

0

z

x

0.5

0.5

11

x

11

(d)

(c)

Fig. 8. Solutions of Problem 3: (a) exact, (b) CDS scheme on uniform grid, (c) HOC scheme on uniform grid, and (d) HOC scheme on nonuniform grid ðk ¼ 0:8Þ.

Table 6 The number of multigrid V-cycles and CPU time with two schemes and different relaxation methods with 323 for Test problem 3. Scheme

k

Point GS

Red–black GS

Four-color GS

Plane GS

Iter.

CPU

Iter.

CPU

Iter.

CPU

Iter.

CPU

CDS

0.0 0.2 0.4 0.6 0.8 0.9

12 12 12 12 18 24

0.547 0.563 0.563 0.562 0.718 0.922

9 11 13 17 28 44

0.469 0.546 0.594 0.718 1.000 1.453

9 12 17 31 – –

0.484 0.594 0.734 1.141 – –

8 9 9 11 12 14

0.531 0.578 0.578 0.641 0.672 0.750

HOC

0.0 0.2 0.4 0.6 0.8 0.9

9 9 10 10 14 18

0.859 0.891 0.953 0.953 1.219 1.531

8 10 11 12 18 21

0.829 0.953 1.031 1.094 1.500 1.687

9 9 10 14 22 24

0.860 0.875 0.937 1.217 1.750 1.875

8 9 10 11 12 14

1.063 1.188 1.313 1.360 1.453 1.656

Y. Ge et al. / Journal of Computational Physics 234 (2013) 199–216

215

grid is also distorted near boundary layers. But the HOC scheme on nonuniform grid with k ¼ 0:8 can get very accurate solution. Table 6 lists the number of multigrid V-cycles with different relaxation methods and the corresponding CPU time in seconds for different stretching parameter k for Test problem 3 with the CDS scheme and the HOC scheme. We find that on uniform grids (k ¼ 0:0) the red–black Gauss–Seidel relaxation takes less CPU time and thus is more cost-effective than the other three relaxations. For all k, the alternating plane Gauss–Seidel relaxation is the most efficient smoother with the least multigrid V-cycle numbers for the two schemes. But in terms of CPU time, for the HOC scheme, when k is small; e.g., k ¼ 0:2; 0:4 (small stretching ratio for grids), the four-color Gauss–Seidel relaxation takes the least time is the most cost-effective. While for big k; e.g., k ¼ 0:6; 0:8; 0:9 (big stretching ratio for grids), the point Gauss–Seidel relaxation is the most cost-effective. For the CDS scheme, when k ¼ 0:0; 0:2, the red–black Gauss–Seidel relaxation takes the least time is the most cost-effective; when k ¼ 0:4; 0:6, the point Gauss–Seidel relaxation is the most cost-effective; and when k ¼ 0:8; 0:9, the alternating plane Gauss–Seidel relaxation is the most cost-effective. 5. Concluding remarks In this paper, a transformation-free HOC finite difference scheme on nonuniform grids for solving the 3D Poisson equation is proposed. The HOC scheme has third- to fourth-order accuracy and produces more accurate solution than the CDS scheme. Numerical studies indicate that for isotropic problems, the HOC scheme, as well as the CDS scheme can achieve its optimal accuracy order on uniform grids. However, for boundary layer problems, fourth-order accuracy for the HOC scheme and second-order accuracy for the CDS scheme are achieved only on nonuniform grids with suitable grid stretching ratios. So, it reveals that grid distributions have significant effects on the computed accuracy of the schemes for boundary layer problems. With reasonable grid distribution on nonuniform grids, we can get more accurate solution than on uniform grids. A multigrid method is employed to solve the linear systems arising from the HOC scheme. The highlight of this multigrid method is that new restriction and interpolation operators on the 3D nonuniform grid are constructed using the volume law. Numerical experiments demonstrate that the present HOC scheme and multigrid method perform well to solve the 3D Poisson equations on nonuniform grids and the computed results show that the present method can dramatically improve the computed accuracy by clustering many more grid points into the boundary layers. Smoothing effects of different relaxation operators are compared and the computed results indicate that for isotropic problems, the pointwise Gauss–Seidel relaxations are good smoothers and the red–black (for the CDS scheme) and the four-color Gauss–Seidel relaxations (for the HOC scheme) are more efficient and cost-effective, but for boundary layer problems, the alternating plane Gauss–Seidel relaxation is more efficient than the other three pointwise Gauss–Seidel relaxations. The present method can also be extended to the design of HOC difference schemes and multigrid method with nonuniform grids for solving other 3D partial differential equations, such as the convection diffusion equations, Helmholtz equations and the Navier–Stokes equations. Acknowledgements This work is supported by the National Natural Science Foundation of China under Grants 11061025 and 11161036, the Key Project of China Ministry of Education under Grant 210239, and the Fok Ying-Tong Education Foundation of China under Grant 121105. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17]

R. Temem, Navier–Stokes Equations, Theory and Numerical Analysis, Elsevier, New York, 1984. C.A.J. Fletcher, Computational Techniques for Fluid Dynamics, Springer, New York, 1991. A.D. Pierce, Acoustics: An Introduction to Its Physical Principles and Applications, Acoustical Society of America, Woodbury, New York, 1989. Y. Jaluria, K.E.Torrance, Computational Heat Transfer, Taylor & Francis, New York, 2003. I. Harari, E. Turkel, Accurate finite difference methods for time-harmonic wave propagation, J. Comput. Phys. 119 (1995) 252–270. U. Schumann, R.A. Sweet, Fast Fourier transforms for direct solution of Poisson’s equation with staggered boundary conditions, J. Comput. Phys. 75 (1988) 123–127. A. McKenney, L. Greengard, A. Mayo, A fast Poisson solver for complex geometries, J. Comput. Phys. 118 (1995) 348–355. U. Schumann, R. Sweet, A direct method for the solution of Poisson equation with Neumann boundary conditions on a staggered grid of arbitrary sizes, J. Comput. Phys. 20 (1976) 171–182. B.L. Buzbee, G.H. Golub, C.W. Nielson, On direct methods for solving Poisson’s equations, SIAM J. Numer. Anal. 7 (1970) 627–656. P.N. Swarztrauber, R.A. Sweet, The Fourier and cyclic reduction methods for solving Poisson’s equation, in: J.A. Schetz, A.E. Fuhs (Eds.), Handbook of Fluid Dynamics and Fluid Machinery, John Wiley and Sons, New York, 1996. T. Chan, R. Glowinski, J. Pe´riaux, O. Widlund, Domain Decomposition Methods, Society for Industrial and Applied Mathematics, Philadelphia, 1989. T. Chan, D.C. Resasco, A domain-decomposed fast Poisson solver on a rectangle, SIAM J. Sci. Stat. Comput. 8 (1987) 27–42. A. Brandt, Multi-level adaptive solution to boundary value problems, Math. Comput. 31 (1977) 333–390. R.A. Nicolaides, On some theoretical and practical aspects of multigrid methods, Math. Comput. 33 (1979) 933–952. M.M. Gupta, J. Kouatchou, J. Zhang, Comparison of second-order and fourth-order discretization for multigrid Poisson solvers, J. Comput. Phys. 132 (1997) 226–232. J. Zhang, Fast and high accuracy multigrid solution of the three dimensional Poisson equation, J. Comput. Phys. 143 (1998). 449-161. J. Zhang, Multigrid method and fourth-order compact scheme for 2D Poisson equation with unequal mesh-size discretization, J. Comput. Phys. 179 (2002) 170–179.

216

Y. Ge et al. / Journal of Computational Physics 234 (2013) 199–216

[18] Y. Ge, Multigrid method and fourth-order compact difference discretization scheme with unequal meshsizes for 3D Poisson equation, J. Comput. Phys. 229 (2010) 6381–6391. [19] M.M. Gupta, J. Kouatchou, J. Zhang, A compact multigrid solver for convection–diffusion equations, J. Comput. Phys. 132 (1997) 123–129. [20] J. Zhang, Accelerated multigrid high accuracy solution of the convection-diffusion equation with high Reynolds number, Numer. Methods Partial Differ. Equat. 13 (1997) 77–92. [21] M.M. Gupta, J. Zhang, High accuracy multigrid solution of the 3D convection-diffusion equation, Appl. Math. Comput. 113 (2000) 249–274. [22] J. Zhang, L. Ge, J. Kouatchou, A two colorable fourth-order compact difference scheme and parallel iterative solution of the 3D convection diffusion equation, Math. Comput. Simulat. 54 (2000) 65–80. [23] D. Bai, A. Brandt, Local mesh refinement multilevel techniques, SIAM J. Sci. Stat. Comput. 8 (1987) 109–134. [24] R. Teigland, I.K. Eliassen, A multiblock/multilevel mesh refinement procedure for CFD computations, Int. J. Numer. Methods Fluids 36 (2001) 519–538. [25] J. Zhang, H. Sun, J.J. Zhao, High order compact scheme with multigrid local mesh Refinement procedure for convection diffusion problems, Comput. Methods Appl. Mech. Eng. 191 (2002) 4661–4674. [26] E. Kalnay De Rivas, On the use of nonuniform grids in finite-difference equations, J. Comput. Phys. 10 (1972) 202–210. [27] L. Farnell, Solution of Poisson equations on a nonuniform grid, J. Comput. Phys. 35 (1980) 408–425. [28] X. Wang, Z. Yang, G. Huang, B. Chen, A high-order compact difference scheme for 2D Laplace and Poisson equations in non-uniform grid systems, Commun. Nonlinear Sci. Numer. Simul. 14 (2009) 379–398. [29] J.C. Kalita, A.K. Dass, D.C. Dalal, A transformation-free HOC scheme for steady convection -diffusion on non-uniform grids, Int. J. Numer. Methods Fluids 44 (2004) 33–53. [30] M.M. Gupta, R.P. Manohar, J.W. Stephenson, High-order difference schemes for two-dimensional elliptic equations, Numer. Methods Partial Differ. Equat. 1 (1985) 71–80. [31] J.Y. Choo, D.H. Schultz, A high order difference method for steady state Navier–Stokes equation, Comput. Math. Appl. 27 (1994) 105–119. [32] W.F. Spotz, G.F. Carey, Formulation and experiments with high-order compact schemes for nonuniform grids, Int. J. Numer. Methods Heat Fluid Flow 8 (1998) 288–297. [33] L. Ge, J. Zhang, High accuracy iterative solution of convection diffusion equation with boundary layers on nonuniform grids, J. Comput. Phys. 171 (2001) 560–578. [34] Y. Ge, F. Cao, Multigrid method based on the transformation-free HOC scheme on nonuniform grids for 2D convection diffusion problems, J. Comput. Phys. 230 (2011) 4051–4070. [35] J. Wang, W. Zhong, J. Zhang, A general meshsize fourth-order compact difference discretization scheme for 3D Poisson equation, Appl. Math. Comput. 183 (2006) 804–812. [36] W. Hackbusch, U. Trottenberg, Multigrid Methods, Springer-Verlag, Berlin, 1982. [37] P. Wesseling, An Introduction to Multigrid Methods, Wiley, Chichester, 1992. [38] C. Liu, Multilevel adaptive methods in computational fluid dynamics, Ph.D thesis, University of Colorado at Denver, 1989. [39] Y. Saad, A flexible inner-outer preconditioned GMRES algorithm, SIAM J. Sci. Stat. Comput. 14 (1993) 461–469. [40] J. Zhang, Preconditioned iterative methods and finite difference schemes for convection–diffusion, Appl. Math. Comput. 109 (2000) 11–30. [41] C.W. Oosterlee, A GMRES-based plane smoother in multigrid to solver 3-D anisotropic fluid flow problems, J. Comput. Phys. 130 (1997) 41–53. [42] Y. Wang, J. Zhang, Fast and robust sixth-order multigrid computation for the three-dimensional convection–diffusion equation, J. Comput. Appl. Math. 234 (2010) 3496–3506.