不可壓縮流的有限元後驗估計與多重網格計算

(1)

行政院國家科學委員會專題研究計畫成果報告

不可壓縮流的有限元後驗估計與多重網格計算

計畫類別：個別型計畫

計畫編號： NSC94-2119-M-009-003-

執行期間： 94 年 11 月 01 日至 95 年 07 月 31 日

執行單位：國立交通大學應用數學系(所)

計畫主持人：吳金典

報告類型：精簡報告

報告附件：出席國際會議研究心得報告及發表論文

處理方式：本計畫可公開查詢

中華民國 95 年 10 月 27 日

(2)

行政院國家科學委員會補助專題研究計畫

v 成果報告 □期中進度報告

（計畫名稱）

不可壓縮流的有限元後驗估計與多重網格計算

計畫類別：v 個別型計畫 □ 整合型計畫

計畫編號：NSC 94－2119－M－009－003－

執行期間：2005 年 11 月 1 日至 2006 年 07 月 31 日

計畫主持人：吳金典

共同主持人：

計畫參與人員：

成果報告類型(依經費核定清單規定繳交)：ｖ精簡報告 □完整

報告

本成果報告包括以下應繳交之附件：

□赴國外出差或研習心得報告一份

□赴大陸地區出差或研習心得報告一份

ｖ出席國際學術會議心得報告及發表之論文各一份

□國際合作研究計畫國外研究報告書一份

執行單位：

國立交通大學應用數學系

中華民國９５年１０月２７日

I

(3)

（一）

計畫中文摘要。（五百字以內）

在本計劃內，我們將發展一個解三維 Navier-Stokes方程有效能和準

確的數值方法來模擬複雜區域內不可壓縮流的行為。為了解

Navier-Stokes 方程, 首先，我们用隱式的 Euler 方法與

Crank-Nicolson 方法來對時間做離散. 因此而產生的非線性對流項

則以 Picard 迭代法或者 Simo-Amero 方法將其線性化. 下一步,

再以 Streamline-Diffusion 的有限元素法 (SDFEM) 離散空間區域.

其中簡單的 P1/P 1元素分別被用來離散速度和壓力兩個變量. 此離

散法由於流線方向有額外的穩定度，使得離散所得的線性系統也因此

比一般由 Galerkin 有限元方法離散所得的線性系統更加的穩定.

再來，我們將集中精力探討怎樣有效地求得這個線性系統的解. 在本

計劃中, 兩種方法將被考慮. 首先,我們考慮一個常用的方法, 在那

裡, 一步的 block Gauss-Seidel 迭代被用來作為線性系統的

preconditioner. 在每一個 block 中則以 GMRES 或 MINRES 來求

解並以多重網格法來改善系統的條件數以增加 GMRES 和 MINRES的

收斂速度. 進而縮短數值計算所需的時間. 另一個方法, 我們分別

于線性系統的左右兩邊乘上特殊的矩陣來改善系統的條件數. 對所

(4)

得的新系統再以 GMRES 迭代法求解. 我們發現第 2 種方法的優勢

在於: (i) 每個 block 中我們只需要一個類似 Poisson 的 solver,

(ii) 計算過程中只要時間的離散尺度和網格的大小固定, 在 (i)

中類似 Poisson 的矩陣並不需隨時間改變而重新計算, (iii) 代數

性的多重網格解是可以使用的, 也就是說我們不需產生粗網格. 因

此程式的複雜度也大大的簡化了.

最後，為了增加數值解的準確度，自適應網格和隨網格大小自動調

整的時間尺度也將一並被考慮. 我們也將用許多標準的流體測試問

題 (包括 driven cavity flow, flow in a backward-facing step 和

flow around a cylinder 等), 來檢測我們的不可壓縮流數值模擬器

的準確性及穩定性.

(二) 計畫英文摘要。（五百字以內）

In this proposal, we would like to develop a robust and accurate

3D Navier-Stokes (NS) solver for simulating the incompressible

flow in a complex domain. For solving the NS equation, first,

both the implicit backward Euler scheme and the 2

nd

order

(5)

Crank-Nicolson scheme will be employed on time domain. The

nonlinear convection term is linearized and treated by Picard

iterations or the Simo-Amero scheme. Next, the

streamline-diffusion finite element method (SDFEM) will be

used for discretization on the spatial domain, where the simple

P1/P1 element is used to discretize the velocity and pressure

variables, respectively. Due to the extra stabilization along

the streamline direction, the resulting bilinear form is

consistent. As a result, the discrete linear system is stable.

Then, we will focus on how to solve the linear system

efficiently. Two approaches will be considered. First, we

consider a heuristic approach where one step of block

Gauss-Seidel iteration is used as a preconditioner of the

system. The multigrid-preconditioned GMRES and MINRES are then

used in the block solvers to accelerate the overall computation

time. Next, our two preconditioners are applied to the linear

system, one on the left and one on the right. The resulting

system is again solved by GMRES. The advantages of the second

(6)

approach includes (i) only block solver which is similar to the

Poisson solver is needed, (ii) no reassembling of matrices is

needed for this block solver during the time evolution, as long

as the time step and the mesh are fixed, and (iii) algebraic

multigrid solver can be applied, i.e. no coarse grids need to

be generated a priori.

Finally, in order to increase the solution accuracy, adaptive

mesh refinement and adaptive time-stepping strategies will be

considered. Many benchmark flow problems including driven

cavity flows, flows in a backward-facing step and flows around

a cylinder will be tested to demonstrate the robustness and

accuracy of the developed solver.

(三) 計畫執行成果:

見下列附件．

(7)

ARTICLE IN PRESS

APNUM:1935

To cite this article: Ming-Chih Lai, et al., An efficient semi-coarsening multigrid method for variable diffusion problems in cylindrical coordinates, Applied Numerical Mathematics (2006), doi:10.1016/j.apnum.2006.07.019

JID:APNUM AID:1935 /FLA [m3SC+; v 1.64; Prn:31/07/2006; 11:30] P.1 (1-10)

Applied Numerical Mathematics••• (••••) •••–•••

www.elsevier.com/locate/apnum

An efficient semi-coarsening multigrid method for variable diffusion

problems in cylindrical coordinates

Ming-Chih Lai

∗

, Chin-Tien Wu, Yu-Hou Tseng

Department of Applied Mathematics, National Chiao Tung University, 1001, Ta Hsueh Road, Hsinchu 300, Taiwan

Abstract

In this paper, we present an efficient multigrid (MG) algorithm for solving the three-dimensional variable coefficient diffusion equation in cylindrical coordinates. The multigrid V-cycle combines a semi-coarsening in azimuthal direction with the red-black Gauss–Seidel plane (radial-axial plane) relaxation. On each plane relaxation, we further semi-coarsen the axial direction with red-black line relaxation in the radial direction. We also prove the convergence of two-level MG with plane Jacobi relaxation. Numerical results show that the present multigrid method indeed is scalable with the mesh size.

Keywords: Multigrid method; V-cycle; Variable diffusion equation; Cylindrical coordinates

1. Introduction

The variable coefficient elliptic equation in three-dimensional domain arises in many physical applications. The heat transfer in heterogeneous material where the thermal conductivity depends on the position is one of the classical examples. The fluid flows with non-constant viscosity is another application. Very often, the geometry we consider is no longer Cartesian, say a circular cylinder instead. Thus, it is more convenient to write the equation in cylindrical coordinates.

In this paper, we consider the following variable coefficient diffusion problem written in cylindrical coordinates on a 3D domain Ω= {(r, θ, z) | 0 < r < 1, 0 θ < 2π, 0 < z < 1} as −1 r ∂ ∂r rβ∂u ∂r + ∂ ∂θ β r ∂u ∂θ + r ∂ ∂z β∂u ∂z = f (r, θ, z) in Ω, (1) u(1, θ, z)= g(θ, z), u(r, θ,0)= p(r, θ), u(r, θ,1)= q(r, θ). (2) Here the diffusion coefficient β(r, θ, z) is inhomogeneous and satisfies 0 < μ β(r, θ, z) 1. The upper bound is not a restriction of the method discussed in this paper since we can re-scale both sides of Eq. (1) by the maximum of β.

When we solve Eq. (1) numerically, the first issue called the coordinate singularity arises. This is because the equation is not valid at r= 0 when it is written in cylindrical coordinates. In [8], Lai et. al. have developed a FFT-based fast direct solver for the Poisson equation (a special case of β(r, θ, z)= 1 in Eq. (1)). The solver relies on the truncated

* _{Corresponding author.}

E-mail address: [email protected] (M.-C. Lai).

(8)

ARTICLE IN PRESS

APNUM:1935

2 M.-C. Lai et al. / Applied Numerical Mathematics••• (••••) •••–•••

Fourier series expansion, where the partial differential equations of Fourier coefficients are solved by the standard centered difference scheme under a shifted grid in the radial direction. The method handles the coordinate singularity without special treatment and the resultant matrix equations can be solved by the fast direct solver such as those provided in public software package—FISHPACK [1]. However, when the diffusion coefficient is inhomogeneous, solving the numerical solution for Eq. (1) is a different story. Since now the elliptic equation has a variable diffusion coefficient, we are unable to write the solution as Fourier series expansion, which means the fast Fourier transform (FFT) cannot be called directly. Furthermore, the equation is not a separable type, the resultant linear system after the discretization cannot be solved by the fast direct solver as mentioned before.

Multigrid (MG) methods are known to be very efficient for solving the elliptic problems. Its main idea consists of applying simple relaxations on the fine grid to smooth the error and correcting the smooth error on the coarse grid, on which the smooth error is relatively oscillatory and can be eliminated more effectively by relaxations. It has been shown that the convergence rates of traditional MG methods are independent of the mesh size for many elliptic problems [2,3]. However, for solving problems with strong anisotropic diffusion or solving problems on anisotropic meshes, the convergence of the MG methods is seriously deteriorated [4,7,13]. The reasons contributing to the slow multigrid convergence in such problems are, first, the traditional Jacobi and Gauss–Seidel relaxations fail to smooth the errors in some directions of the domain and, second, the standard coarse grids cannot represent the errors that are highly frequent in these directions. Many cures have been proposed to improve the convergence rate. For some two-dimensional anisotropic problems on the Cartesian coordinates, Brandt, Hackbush and Mulder [4,7,11] have obtained robust convergence of multigrid in which coarse grids from semi-coarsening are employed. For three-dimensional anisotropic problems, Dendy [6], Llorente and Melson [10], and Washio and Oosterlee [12] obtained robust conver-gence by using semi-coarsened grids and the alternating plane relaxations in their MG methods. Schaffer [14] also developed an efficient semi-coarsening multigrid method for symmetric and nonsymmetric elliptic PDEs with highly discontinuous and anisotropic coefficients in two- and three-dimensional Cartesian domains. In this paper, we extend the discretization used in [9] to the 3D variable diffusion equation in cylindrical coordinates (1)–(2), and present an efficient semi-coarsening MG algorithm to solve the resultant linear system.

The rest of this paper is as follows. In Section 2, we introduce the finite difference discretization to Eq. (1). In Section 3, we first review some elements of the multigrid method and then provide a two-level convergence analysis for the semi-coarsening MG with Jacobi plane relaxation. The numerical results including the accuracy check and the detailed performance comparison for the present multigrid algorithm are shown in Section 4.

2. Finite difference discretization

We define the same grid points in the radial, azimuthal and axial directions as in [8] by

ri = (i − 1/2)r, ri+1/2= ri+ r/2, ri−1/2= ri− r/2, (3)

θj= jθ, θj+1/2= θj+ θ/2, θj−1/2= θj− θ/2, (4)

zk= kz, zk+1/2= zk+ z/2, zk−1/2= zk− z/2, (5)

where r= 2/(2nr+ 1), θ = 2π/nθ and z= 1/(nz+ 1). By the choice of the radial mesh width, the boundary

values are defined on the grid points. Let the discrete values be denoted by uij k ≈ u(ri, θj, zk), fij k≈ f (ri, θj, zk),

gj k≈ g(θj, zk), pij ≈ p(ri, θj)and qij ≈ q(ri, θj). The values at the half-integered points are defined in a similar

fashion, such as βi+1/2,j,k≈ β(ri+1/2, θj, zk). Using the centered difference method to discretize Eq. (1), we have

− ri+1/2βi+1/2,j,k ui+1,j,k− ui,j,k r − ri−1/2βi−1/2,j,k ui,j,k− ui−1,j,k r r + βi,j_+1/2,k ri ui,j_+1,k− ui,j,k θ − βi,j_−1/2,k ri ui,j,k− ui,j_−1,k θ θ + ri βi,j,k_+1/2 ui,j,k+1− ui,j,k z − βi,j,k−1/2 ui,j,k− ui,j,k−1 z z = rifi,j,k. (6)

Among the above representation, the numerical boundary values are given by unr+1,j,k = gj,k, ui,j,0 = pi,j, ui,j,nz+1= qi,j, and ui,0,k= ui,nθ,k, ui,nθ+1,k= ui,1,k since u is 2π periodic in the azimuthal direction. At i= 1,

(9)

ARTICLE IN PRESS

APNUM:1935

JID:APNUM AID:1935 /FLA [m3SC+; v 1.64; Prn:31/07/2006; 11:30] P.3 (1-10) M.-C. Lai et al. / Applied Numerical Mathematics••• (••••) •••–••• 3

we have immediately observed from Eq. (3) that r1/2= 0, so the coefficient of u0,j,k is zero. This implies that the

scheme does not need any extrapolation for the inner numerical boundary value u0,j,kso that there is no pole condition

needed.

Let us order the unknowns to a solution vector u, the resultant linear system of (6) can be written as Au= f . Here the matrix A is given by

A= Iθ⊗ Brz+ Cθ⊗ Irz, (7)

where⊗ is the regular Kronecker product, Brzis the matrix with a five-point stencil obtained from discretization of

the operator _∂r∂ (rβ∂u_∂r)+_∂z∂(rβ∂u_∂z), Cθ is the circulant matrix with a three-point stencil obtained from discretization

of the operator _∂θ∂(β_r ∂u_∂θ), and Iθ and Irz are the identity matrix of size nθ × nθ and nr · nz× nr · nz, respectively.

Clearly, in the region close to the origin, the problem (1) is strongly anisotropic, namely, the diffusion in the azimuthal direction is much greater than that in other directions. As a result, the matrix A becomes indefinite as the grid space decreases. Therefore, a scalable and efficient solver for the linear system Au= f is desired.

3. Semi-coarsening multigrid method

In this paper, we would like to solve the aforementioned linear system by the semi-coarsening MG with plane relaxations. In the following, firstly, the general multigrid V-cycle algorithm is introduced, and then the convergence of a two-level V-cycle of the MG method is proved under the assumption that the solution of Eq. (1) is smooth enough so that the discretization error of the scheme (6) is second-order accurate.

Given a sequence of meshes{G}, let Vbe the vector space of the nodal variables on G, Athe matrix obtained from the difference equation (6) on G, and wthe initial guess. Let MG1represents the direct solver on the coarsest

grid G1. A typical multigrid V-cycle on -level is shown recursively in Algorithm 1. Here, fis the right-hand side obtained from (6) on G, (M)−1 represents the smoothing operator and the operators I−1(restriction) and I₋₁ (prolongation) represent the grid transfers between G−1_{and G}_{. The elements of the multigrid V-cycle are described}

as follows.

Algorithm 1. Multigrid V-cycle MG(w, f)

(1) Set u= w;

(2) (Pre-smoothing) u= u+ (M)−1(f− Au), ν1times;

(3) (Restriction) ¯f= I−1(f− Au);

(4) (Coarser grid correction) q1= MG−1(0, ¯f);

(5) (Prolongation) ¯q1= I−1q1, and set u= u+ ¯q1;

(6) (Post-smoothing) u= u+ (M)−1(f− Au) ν2, times;

(7) Set w= MG(w, f)= u.

For the MG considered here, (i) the coarse grids G−1are obtained from semi-coarsening the fine grid G along the azimuthal direction, (ii) the prolongation operator I₋₁is the linear interpolation along the azimuthal direction and the restriction operator I−1is the transpose of the interpolation, (iii) the matrix A−1is obtained from discretization scheme (6) on (− 1)-level grid G−1and (iv) the smoothing operator M consists of the Jacobi, or Gauss–Seidel plane relaxations. From now on, we use the notations MGJ(ν1, ν2)and MGGS(ν1, ν2)to represent the different

ver-sions of the above multigrid method in which ν1times of pre-smoothing and ν2times of post-smoothing by plane

Jacobi and plane Gauss–Seidel smoothers, respectively, are applied.

3.1. The convergence of two-level MG

In the following, we present the convergence analysis for the two-level V-cycle MG method with Jacobi plane relaxations. To simplify our convergence proof, we assume the number of grid points used in r, θ and z directions are

nr= nθ/2= nz= n. Let · 2 denote the standard vector 2-norm or matrix 2-norm depending on which quantities

are considered. LetwQ=

√

(10)

ARTICLE IN PRESS

APNUM:1935

positive definite matrix and (·, ·) is the inner product on V_{. For any linear transformation T from V}_{to V}_{, the matrix}

Q-norm, denoted byT Q, is the matrix 2-norm associated with the Q-weighted Euclidean norm.

First, let us show that the Jacobi plane relaxation is convergent. Recall that the Jacobi plane relaxation is obtained from the matrix splitting

A= M− N, (8)

where

M= Iθ⊗ Brz+ Dθ⊗ Irz and N= (Dθ− Cθ)⊗ Irz, (9)

and the matrix Dθis the diagonal part of Cθ. The matrices Mand Ncan also be written in the following block form

M= ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ M1 0 0 M2 0 . ._. . ._. . ._. 0 M2n−1 0 0 M2n ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , N= ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 N1 0 0 N2n N1 0 N2 0 0 0 N2 . .. . .. ... 0 0 . .. 0 N2n₋₁ N2n 0 0 N2n−1 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , (10)

where Mj are the block tridiagonal matrices with the stencils

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ −ri−1/2βi−1/2,j,k r2 , ri−1/2βi−1/2,j,k+ri+1/2βi+1/2,j,k r2 + riβi,j,k−1/2+riβi,j,k+1/2 z2 + βi,j−1/2,k+βi,j+1/2,k riθ2 , −ri+1/2βi+1/2,j,k r2

i,k in the main diagonal blocks,

−riβi,j,k−1/2

z2

i,k in the upper off-diagonal blocks,

−riβi,j,k+1/2

z2

i,k in the lower off-diagonal blocks,

(11)

and the matrices Nj, j = 1 . . . 2n, are the block diagonal matrices with stencils [ βi,j+1/2,k

riθ2 ]i,k, respectively. Let S

J =

(M)−1N be the error reduction operator of the Jacobi plane relaxation. It is clear that, due to the given Dirichlet boundary conditions, Mand Aare irreducible M-matrices [15]. Therefore, (M)−1and (A)−1are positive. Since

N is non-negative, the splitting (8) is a regular splitting. As a result, the Jacobi plane relaxation SJ is convergent

sinceS_J2<1 ([15] Theorem 3.29).

Throughout this section, we define the matrix Ras

R= Iθ⊗ Rrz = Iθ⊗

Iz⊗ Rr

, (12)

where R_r= diag(ri). Obviously, one also hasSJR<1.

Lemma 3.1. Suppose Mj, j= 1, . . . , 2n and Nl, l= 1, . . . , 2n are the block matrices defined in (10). The inequality

M_j−1Nl_R rz<

κ

2μ (13)

holds, where κ < 1 is a constant depending on β(r, θ, z) and μ is the lower bound of β(r, θ, z).

Proof. The matrix Mjcan be rewritten as

Mj=

1

θ2

R_rz−1Bj(Irz+ C),

where Bj is a diagonal matrix whose diagonal entries [ ¯βj]i,k are defined by ¯βj = βi,j−1/2,k+ βi,j+1/2,k, and the

matrix C= θ2B_j−1R_rz Mj− Irz. Clearly, one has

M_j−1Nl_R rz Irz+ R_rz1/2CR_rz−1/2−1θ2B_j−1Nl₂ (14) =Irz+ R_rz1/2CR_rz−1/2−1 βi,l+1/2,k ¯βj i,k 2 . (15)

(11)

ARTICLE IN PRESS

APNUM:1935

From stencils (11) and the given Dirichlet boundary conditions in (2), one can conclude that C is an irreducibly di-agonal dominant matrix. Therefore, the matrix C is positive definitive. Moreover, (R_rz)1/2C(R_rz)−1/2is also positive definite. Let Λmin(C)denote the smallest eigenvalue of the matrix C= (Rrz )1/2C(Rrz)−1/2. Obviously, the

assump-tion that μ < β 1 and (14) imply M_j−1Nl_R rz 1 1+ Λmin(C) βi,l+1/2,k ¯βj i,k 2 (16) < κ 2μ. (17)

Hence, the inequality (13) holds. 2

Corollary 3.2. For the constant diffusion case, suppose the mesh width satisfies the relation as r= η1z= η2θ 1, where η1=2n_2nz+2

r+1 and η2=

nθ

(2nr+1)π are the grid aspect ratios of z and θ with respect to r, and η1≈ η2≈ 1. Then, we have

M_j−1Nl_R rz

κ

2, (18)

and κ≈ (1 − cr2) for some constant c >0 independent of mesh size. As a result of (18), the plane Jacobi relaxation

converges with

S_J

R 1 − cr2. (19)

Proof. First, the stencil of the matrix C in the proof of Lemma 3.1 can be rewritten as C= θ r 2 C, (20)

where Cis a block tridiagonal matrix with[−ri√riri−1/2, (2+ 2η2₁)r_i2,−η2₁ri√riri+1/2] in the diagonal blocks, and

[−r2

i] and [−η21ri2] in the lower and upper off-diagonal blocks, respectively. By the Gersgorin theorem, we obtain the

following estimates as Λmin(C)− 2+ 2η₁2r_i2 ri 1+ η2₁ri+ √riri−1/2+ η21 √_r iri+1/2 2+ 2η21

r_i2− c0rir, for some constant c0>0

2+ 2η2₁r_i2− 0.5c0r2. (21)

Clearly, from (20) and (21), there exists a constant c1≈ 0.5c0(θ_r)2= 0.5c0(_η1₂)2>0 such that Λmin(C) c1r2.

Next, let e1= [1, 0, . . . , 0]. Since

Ce12= θ r 2 (2+ 2η1)r12 2 + η2 1r 2 2r2r3/2+ r14 1/2 = 16η2₁+ 60η1+ 5 4η₂2 r 2_{= c} 2r2

from the positive definiteness of C, one can conclude that Λmin(C) c2r2, for some constant c2>0. By using the

upper and lower bounds of Λmin(C), the inequality (16) can now be shown easily as

M_j−1Nl_R rz 1 2(1+ Λmin(C)) 1 2 1− Λmin(C)+ Λmin(C)2 1 2 1− c1r2+ c22r4 1 2 1− cr2,

for some constant c > 0. The last inequality arises easily from the assumption that the radial mesh r is much smaller than one. The inequality (19) can now be derived directly from the Gersgorin theorem. 2

(12)

ARTICLE IN PRESS

APNUM:1935

Now, we are ready to prove the convergence of the 2-level MG. First, the error reduction operator of the 2-level MG for solving the linear system obtained from (6) can be written as

EMG=

A−1− I₋₁A−1−1I−1RAS_Jν,

where ν is the number of relaxations. To prove the convergence of MG is independent of mesh size, we need to show that there is some number ν > 0 such that

EMGR<1. (22)

In the following, we let e _{be an arbitrary error vector in V} _satisfying_e

R= 1, and (f)s = A(S_J)νe be the

corresponding residual of the relaxed error. Clearly, one has

RA−1− I₋₁(A−1)−1I−1RA(S_J)νe, e

=RR−1A−1− I₋₁R−1−1A−1−1I−1fs, e

R−1A−1− I₋₁R−1−1A−1−1I−1fs

R.

Since (R)−1Aand (R−1)−1A−1represent the discrete difference operators in (6) on grids G and G−1, respec-tively, the truncation errors are O(r2)(the meshes θ and z are the same order as r) if the solution is smooth. Thus, there is a constant c3>0 such that

EMGR c3r2sup e RAS_Jνe, AS_Jνe1/2 c3r2A S_Jν R c3r2ASJRSJ ν₋₁ R. (23)

Next, since AS_J= N(I− S_J)and r2N_R 2 βi,l+1/2,k π2 i,k 2 2 π2,

clearly, the inequality (23) implies EMGR 2c3 π2I− S JRS J ν−1 R 4c3 π2S J ν−1 R . (24)

Therefore, from the factS_J_R<1, the inequality (22) holds when the number of Jacobi plane relaxations is large

enough. We conclude our result of the multigrid convergence in the following theorem.

Theorem 3.3. For the variable diffusion problem (1) discretized by the finite difference scheme (6) on the cylindrical

uniform grids G, let S_J and EMGbe the error reduction operator of plane Jacobi relaxation and the present 2-level MG, respectively. There exists a positive constant c4independent of the grid size such that

EMGR< c4SJ ν−1

R , (25)

where ν is the number of smoothing steps in the MG. If ν is large enough, the 2-level MG is convergent independent of the mesh size.

Remark 3.4. From the convergence rate of the plane Jacobi relaxation estimated by the inequality (19) of

Corol-lary 3.2, Theorem 3.3 suggests that for the constant diffusion case, the number of Jacobi plane relaxations should increase four times when the mesh is refined to achieve grid-independent convergence of MG. On the other hand, since the plane Jacobi matrix S_J >0 with zero diagonal blocks and (R)1/2S_J(R)−1/2 is block-symmetric, irre-ducible and convergent, by the Perron–Frobenius theorem, and Theorem 4.15 of the Varga’s book [15], the relation of the asymptotic convergence rates between the plane Jacobi S_J and plane Gauss–Seidel S_GS can be derived as follows.

S_GS _R S JR 2− S_JR 1− cr2 1+ cr2≈ 1 − 2cr 2_.

(13)

ARTICLE IN PRESS

APNUM:1935

Fig. 1. The 2-norm of relative residuals with grid sizes n= 16(∗), 32(), 64(2) for MGJ(left) and MGGS(right).

Moreover, the estimate of the error reduction rate of MG with S_√ _GS smoother in (24) becomes EMGR 10c

π2 SGS

ν−1

R . Hence, the grid-independent convergence of MG should be easier to achieve when plane Gauss–

Seidel relaxations are employed in the pre-smoothing and post-smoothing steps of the MG.

Fig. 1 shows the convergence behaviors of MG for the constant diffusion case β= 1. One can see that when the plane Jacobi smoother is employed in MG, in order to have grid-independent convergence rate, the number of smoothing steps must increase four times as the grid number n doubles. This is exactly what we expect from the previous remark. In practice, only one step of the plane Gauss–Seidel relaxation is required for convergence and the convergence rate is independent of the problem size as shown in Fig. 1.

4. Numerical results

In this section, we report several numerical tests for the multigrid method. On a given nr× nz× nθmesh size with

nr= nz= nθ/2, the present multigrid algorithm is a complete log2nθ-level V-cycle which employs semi-coarsening

in the azimuthal (θ ) direction, and one step of red-black plane (r–z plane) Gauss–Seidel relaxation for both the pre-smoothing and the post-pre-smoothing. It is important to note that the present MG method semi-coarsens the grid in the

θ direction instead of the r direction. This is because if the radial grid is coarsened, then the coarser grid would not coincide with the fine grid based on our radial grid arrangement (3). The other reason is that Eq. (1) is on r, θ, z rectangular domain, where the strongly anisotropic diffusion occurs along the θ direction near the z-axis. It is well known that the semi-coarsening direction is determined by which direction is strongly anisotropic [5].

In order to save the computation time, each plane relaxation step is further approximated by one step of log2nz

-level multigrid V-cycle iteration where coarse grids are obtained by semi-coarsening in the axial (z) direction on each plane and employing a red-black line Gauss–Seidel relaxation in the radial (r) direction. The restriction (fine to coarse)

(14)

ARTICLE IN PRESS

APNUM:1935

Fig. 2. The 2-norm of relative residuals of MGRB(1, 1) and the maximum norm of errors for Example 1 (smooth diffusion coefficient) with grid

sizes n= 16(∗), 32(), 64(2), 128(◦).

and prolongation (coarse to fine) operators are the conventional full weighting and linear interpolation, respectively. Hereafter, we denote this version of the multigrid method as MGRB(1, 1).

In this paper, we have tested the convergence behavior and the numerical accuracy of the present MG algorithm by the following examples. In those examples, we choose the analytical solution u(x, y, z)= sin(πx) sin(πy) sin(πz) and vary the diffusion coefficients based on different smoothing characteristics as

(1) β= 1 + sin2(π(x+ y + z)) (smooth case);

(2) β= 1 + sin2(10θ ) (moderately oscillatory in θ direction);

(3) β= 1 + sin2(10θ ) sin2(20π z) (moderately oscillatory in θ and z directions). The right-hand side functions are obtained by substituting the solutions into Eq. (1).

The stopping tolerance of the iterations is set to be 10−8for all of our test cases. All numerical runs are carried out on a PC with 1 GB RAM and 1.3 GHz CPU, and the program is written in MATLAB. All the following figures show the 2-norm of the relative residuals and the maximum norm of the errors versus the number of V-cycles. From Fig. 2, we can see that the residual norms decrease as the same rate for all mesh sizes and it takes about 8 V-cycles to reach the stopping criterion. However, it takes only 5 V-cycles to reach the discretization error for all meshes. Furthermore, the maximum norm of errors decrease by a factor of four as the resolution doubles which indicates that our finite difference discretization for Eq. (1) is indeed second-order accurate. Fig. 3 shows the convergence and accuracy performance for the case of moderately oscillatory diffusion in θ direction. This example takes about 3–4 more V-cycles than the smooth case to reach the residual stopping criterion and the discretization error. It leads to the same convergence rate for all meshes as well. Note that, since the diffusion coefficient is moderately oscillatory, the mesh n= 16 is under-resolving for this case. Fig. 4 shows the performance of our multigrid method for the case of moderately oscillatory diffusion in θ and z directions. Although the meshes n= 16, 32 show under-resolution accuracy, the residual convergence rates are still the same for all meshes. All three examples show our present multigrid V-cycle is scalable with the mesh size perfectly.

(15)

ARTICLE IN PRESS

APNUM:1935

Fig. 3. The 2-norm of relative residuals of MGRB(1, 1) and the maximum norm of errors for Example 2 (diffusion coefficient moderately oscillatory

in θ direction) with grid sizes n= 16(∗), 32(), 64(2), 128(◦).

Fig. 4. The 2-norm of relative residuals of MGRB(1, 1) and the maximum norm of errors for Example 3 (diffusion coefficient moderately oscillatory

(16)

ARTICLE IN PRESS

APNUM:1935

Acknowledgements

The first author was supported by National Science Council of Taiwan under research grant NSC-94-2115-M-009-004. The second author was supported by National Science Council of Taiwan under research grant NSC-94-2119-M-009-003.

References

[1] J. Adams, P. Swarztrauber, R. Sweet, Fishpack—a package of Fortran subprograms for the solution of separable elliptic partial differential equations, 1980. Available in: http://www.netlib.org/fishpack.

[2] J.H. Bramble, J.E. Pasciak, J. Xu, The analysis of multigrid algorithms for nonsymmetric and indefinite elliptic problems, Math. Comp. 51 (1988) 398–414.

[3] J.H. Bramble, D.A. Kwak, J.E. Pasciak, Uniform convergence of multigrid V-cycle iterations for indefinite and nonsymmetric problems, SIAM J. Numer. Anal. 31 (1994) 1746–1763.

[4] A. Brandt, Multi-level adaptive solutions to boundary value problems, Math. Comp. 31 (1977) 333–390. [5] W.L. Briggs, V.E. Henson, S.F. McCormick, A Multigrid Tutorial, second ed., SIAM, Philadelphia, PA, 1999.

[6] J.E. Dendy Jr., Two multigrid methods for three dimensional problems with discontinuous and anisotropic coefficients, SIAM J. Sci. Stat. Comput. 8 (1987) 673–685.

[7] W. Hackbush, Multigrid Methods and Applications, Springer, Berlin, 1985.

[8] M.-C. Lai, W.-W. Lin, W. Wang, A fast spectral/difference method without pole conditions for Poisson-type equations in cylindrical and spherical geometries, IMA J. Numer. Anal. 22 (2002) 537–548.

[9] M.-C. Lai, Y.-H. Tseng, A fast iterative solver for the variable coefficient diffusion equation on a disk, J. Comput. Phys. 208 (2005) 196–205. [10] I.M. Llorente, N.D. Melson, Robust multigrid smoothers for three-dimensional elliptic equations with strong anisotropies,

NASA/CR-1998-208700 ICASE Report No. 98–37.

[11] W.A. Mulder, A high-resolution Euler solver based on multigrid, semi-coarsening, and defect-correction, J. Comput. Phys. 100 (1992) 209– 228.

[12] T. Washio, C. Oosterlee, Flexible multiple semicoarsening for three dimensional singularly pertubered problems, SIAM J. Sci. Comput. 19 (1998) 1646–1666.

[13] A. Reusken, On a robust multigrid solver, Computing 56 (1996) 303–322.

[14] S. Schaffer, A semicoarsening multigrid method for elliptic partial differential equations with highly discontinuous and anisotropic coeffi-cients, SIAM J. Sci. Comput. 20 (1998) 228–242.

(17)

(四) 計畫成果自評：

本計劃已完成二維不可壓縮流的數值計算．其中時間區域已達二階精

度而空間的 Streamline-Diffusion 離散法也使我們能計算高雷諾

數的流體．我們也成功的將多重網格法及自適性網格運用於二維不

可壓縮流的計算．其中多重網格於三維變係數擴散方程的應用已發

表於國際期刊. 其他部分成果請見下三頁．雖然本計劃最終目標雖

尚未完成，但考慮三維的流體計算是非常困難與重要，我們寄望未來

能繼續在國科會的支持下完成未竟之業．

1

(18)

(a) Adaptive meshes for various benchmark problems:

Driven Cavity Flow (Re=10

4

)

Flow around Cylinder (Re=200)

Flow around airfoil NACA0012 at various times

(Re=10

6

, attack angle =20

o

)

(19)

(b) Accuracy are achieved for the benchmark problem (flow

around cylinder)

0.19 0.19(4)(4) 0.2 0.2 0.167 0.167(1)(1) 0.164 0.164(4)(4) 0.167 0.167 0.12~0.13 0.12~0.13(1)(1) 0.139 0.139(3)(3) 0.125 0.125 -ST ST ± ±0.690.69(4)(4) ± ±0.70.7 ± ±0.3390.339(4)(4) ± ±0.350.35 -C C_L_L 1.16 1.16(1)(1) 1.31 1.31±±0.0490.049(4)(4) 1.38 1.38±±0.050.05 1.24 1.24(1)(1) 1.35 1.35±±0.0120.012(4)(4) 1.38 1.38±±0.010.01 1.41 1.41(1)(1) 1.38 1.38(3)(3) 1.46 1.46 2.22 2.22(1)(1) 2.19 2.19(2)(2) 2.11 2.11 C C_D_D 200 200 100 100 50 50 20 20 Re Re

(1) for Re=100 and Re=200 are obtained experimentally by Clift and Reshko

respectively.

(2) is computed on a 640x320 grid by Donna Calhoun, Courant Institute of

mathematics science, in 2002.

(3) is computed on a 267x147 grid by Saki and Biringen, Dept. of Areospace

Engineering, Univ. of Colorado, in 1996.

(4) is computed on a 256x256 grid by Liu, Zheng and Sung in 1998.

In the following, C_D=drag coefficient= C_L= lift coefficient=

ST= Strouhal number= 2

2 F

d

U L

ρ

∞ 2 2F_l U L

ρ

∞ freq L

ν

∗

The numbers in blue color are obtained from our simulations.

Clearly, we have achieved the same accuracy as other

researchers. However, the mesh points generated in our

calculation are just about 7000 to 8000 points which is much

less than the mesh points used in all the other calculation.

(20)

(c) Efficiency is observed for the benchmark problem (flow

around cylinder)

AMG is efficient

Tables on the right show the average number of Iterations for various flows.

200 time steps are performed on each grid. The initial time step is 0.2 and step size is halved when mesh is refined. 21.6 21.6 28.1 28.1 19.3 19.3 GMRES

GMRES--AMGAMG

40.1 40.1 50.5 50.5 19.3 19.3 GMRES GMRES--GMGGMG 4664 4664 2662 2662 1448 1448 nodes nodes Re=200 17.1 17.1 22.4 22.4 16.3 16.3 GMRES

GMRES--AMGAMG

25.6 25.6 29.4 29.4 13.1 13.1 GMRES GMRES--GMGGMG 4664 4664 2662 2662 1448 1448 nodes nodes Re=20

不可壓縮流的有限元後驗估計與多重網格計算

行政院國家科學委員會專題研究計畫 成果報告