Direct Solutions of Finite-Difference Systems for Poisson's Equation II. Complex Cases

Science and

Technology

Copyright © 2013 by JSME

Direct Solutions of Finite-Difference Systems for

Poisson’s Equation II. Complex Cases

∗∗_{Institute of Fluid Science, Tohoku University} Katahira 2–1–1, Aoba-ku, Sendai, Japan

E-mail: [email protected]

Abstract

The direct method developed in the first paper (J. Comp. Sci. Tech. 6 (2012), 147-156) to solve the finite-difference systems for Poisson’s equation is extended to the cases with more complex boundary conditions. The first case is the plasma reactor of the international standard and the second is the three-dimensional plasma reactor. In the first case the computer execution time for the present method is 1 to 3% of the time for SOR with Chebyshev acceleration. In the second case the computer time for the present method is nearly the same as that for FFT method.

Key words : Linear Systems, Band Matrix, Direct Method, Poisson Solver,

Finite-Difference Method

1. Introduction

Electronics industry is requesting PIC/MC simulations(1),(2)_{of materials processing} plas-mas with realistic dimension and complexity.(3)_{The electric field in plasmas is governed by} Poisson’s equation. To answer the request we need a quick and accurate solver of Poisson’s equation subject to complex boundary conditions.

In the previous paper(4)_{we considered the simplest boundary condition for the} axisym-metrical Poisson’s equation. Here we treat more complex cases. We first consider the repre-sentative plasma reactor named GEC cell (Gaseous Electronics Conference cell)(5), that has been widely used among theoreticians and experimentalists as the international standard. Here the finite-difference systems for the axisymmetrical Poisson’s equation for this plasma reactor are solved by extending the previous direct method.(4)_{The computer execution time for the} obtained method is compared with that for SOR with Chebyshev acceleration. Adopting the double precision computation, we used a single processor (class s, 100 GFLOPS) of the vec-tor supercomputer NEC SX-9 in Cyber Science Center of Tohoku University. The processor has a function of automatic vector operation. We employed the function but did no efforts to enhance the vector operation ratio. Two sample calculations were performed : (1) 143000 grid points and (2) 358000 grid points. In case (1) the computer time 0.22s for the present method is 2.6% of the time for SOR and in case (2) the time for the present method is 1.1% of that for SOR.

Secondly, the finite-difference system for three dimensional Poisson’s equation is solved by extending the previous direct method.(4)_{The computational domain is the box with 50×} 50× 50 grids. The finite-difference system is solved by the present method and also by the FFT(fast Fourier-transform) method. The computer time 0.37s for the present method was nearly equal to the time 0.34s for the FFT. Note, however, that the FFT method is applicable only to very simple boundary conditions.

The direct method in the first paper can be extended to problems with geometry more complex than in this paper. Since such a task of making the algorithm may take time, however, it appears to be better to use the algorithm in Section 2 in the first paper, which is applicable in principle to any band matrix resulting from complex geometry.

*

Kenichi NANBU

*Received 11 July, 2013 (No. 13-00138) [DOI: 10.1299/jcst.7.426]

2. Axisymmetrical GEC cell

Figure.1 shows GEC cell. This is a typical plasma reactor for etching or sputtering used in semiconductor industry. Two electrodes with radii R1and R2 are in the circular cylinder with radius R. The surface of the lower electrode is at z= H1and that of the upper electrode is at z= H2. The height of the reactor is H.

In the axisymmetrical coordinates system Poisson’s equation is ∂2_Φ ∂r2 + 1 r ∂Φ ∂r + ∂2_Φ ∂z2 = g(r, z), (1)

where the domain is 0≤ r ≤ R and 0 ≤ z ≤ H. We choose the boundary conditions as follows ∂Φ ∂r =0 at r= 0, (H1≤ z ≤ H2), (2a) Φ = 0 at r = R, (0 ≤ z ≤ H), (2b) Φ = 0 at z = 0, (R1 ≤ r ≤ R), (2c) Φ = 0 at z = H, (R2≤ r ≤ R), (2d) Φ = V1(r) at z= H1, (0 ≤ r ≤ R1), (2e) Φ = V2(r) at z= H2, (0 ≤ r ≤ R2), (2f) Φ = W1(z) at r= R1, (0 ≤ z ≤ H1), (2g) Φ = W2(z) at r= R2, (H2≤ z ≤ H). (2h) If we think ofΦI,JasΦ (IΔr, JΔz), the difference equation takes the form

 1− 1 2I  ΦI−1,J− 2 (1 + β) ΦI,J+  1+ 1 2I  ΦI+1,J + β(ΦI_,J−1+ ΦI_,J+1)= (Δr)2gI,J, (3) 1 2 5 3 4 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84

O

R

R

r

H

H

H

z

R

where β= (Δr/Δz)2and gI,J= g (IΔr, JΔz) .

We emphasize that the FFT method is not applicable in solving (3) because of complexity of the boundary conditions. A sample of grid points and serial grid numbers is shown in Fig.1. We divide the range [0, R] into M+ 1 intervals and the range [0, H] into U + 1 intervals. The grid size isΔr = R/ (M + 1) and Δz = H/ (U + 1). The grid point (I, J) denotes r = IΔr and z= JΔz, where I = 0, 1, ..., M + 1 and J = 0, 1, ..., U + 1. Let us change (I, J) into the serial grid number j. To do this we divide the domain into three regions.

( I ) Lower region: 0≤ z ≤ H1and R1 ≤ r ≤ R. We put L grid points in the r-direction and U1points in the z-direction. The number of grid points is N1(= LU1). (The grid points on the boundaries are not included.)

( II ) Middle region: H1 < z < H2 and 0 ≤ r ≤ R. We put M grid points and U2 grid points in the r- and z-direction respectively. The number of grid points in this region is N2(= MU2).

(III) Upper region: H2≤ z ≤ H and R2≤ r ≤ R. We put N grid points and U3grid points in the r- and z-direction respectively. The number of grid points in this region is N3(= NU3). The total grid points in the domain amount to N1+ N2+ N3(= n). This is the number of unknowns. In the lower region the serial grid number j for the grid point (I, J) is

j= M + 1 − I + L(J − 1),

where j= 1, 2, ..., N1for I = M, M − 1, ..., M − L + 1 and J = 1, 2, ..., U1. Note that j increases with decreasing I. See Fig. 1. In the middle region the serial grid number is

j= N1+ (M + 1 − I) + M(J − U1− 1),

where j= N1+1, N1+2, ..., N1+ N2for I= M, M −1, ..., 1 and J = U1+1, U1+2, ..., U1+U2. In the upper region the serial grid number is

j= N1+ N2+ M + 1 − I + N(J − U1− U2− 1),

where j= N1+ N2+ 1, N1+ N2+ 2, ..., n for I = M, M − 1, ..., M − N + 1 and J = U1+ U2+ 1, U1+ U2+ 2, ..., U.

Let us write xj= ΦI,J. Then we have xj+1= ΦI−1,Jand xj−1= ΦI+1,J. In the lower region (3) takes a general form

bL, jxj−L+ b1, jxj−1+ djxj+ a1, jxj+1+ aL, jxj+L= fj, (4) where fj= (Δr)2gI,Jand

dj= −2(1 + β), (5a) a1, j= 1 − 1 2I, (5b) b1, j= 1 + 1 2I, (5c) aL, j= bL, j= β. (5d)

Some terms in (4) must be removed owing to the boundary conditions of (2b), (2c), and (2g). After the replacement

fj=⇒ fj− a1, jW1(JΔz) , (I = M − L + 1; J = 1, 2, ..., U1),

we remove a1, jxj₊₁ from (4). Equation (2b) requires removal of b1, jxj₋₁ for I = M and

J = 1, 2, ..., U1, and (2c) requires removal of bL, jxj_−L for I = M, M − 1, ..., M − L + 1 and

J= 1. After all (4) becomes

b1, jxj−1+ djxj+ a1, jxj+1+ aL, jxj+L= fj, ( j = 2, 3, ..., L), (6b)

bL, jxj_−L+ b1, jxj₋₁+ djxj+ a1, jxj₊₁

+ aL, jxj+L= fj, ( j = L + 1, L + 2, ..., N1). (6c)

Regard fjas the quantity after correction of the boundary condition. For convenience’ sake we left null a1, jxj₊₁for j= L, 2L, ..., U1L and null b1, jxj₋₁for j= L+1, 2L+1, ..., (U1−1)L+1 in (6). The grid points j+ L for j = N1− L + 1, N1− L + 2, .., N1lie in the middle region.

In the middle region (3) takes a general form

bM, jxj−M+ b1, jxj−1+ djxj+ a1, jxj+1+ aM, jxj+M= fj, (7)

where dj, a1, j, and b1, jare given by (5) and aM, j= bM, j= β. Correction of (7) due to boundary conditions (2a), (2b), (2e), and (2f) is as follows. Equation (2a) requires the replacement

a1, jxj+1=⇒ a1, j  4 3xj− 1 3xj−1  , (I = 1; J = U1+ 1, U1+ 2, ..., U1+ U2). That is, dj=⇒ dj+ 4 3a1, j, b1, j=⇒ b1, j− 1 3a1, j. Other replacement is fj=⇒ fj− bM, jV1(IΔr) , (I = 1, 2, ..., M − L; J = U1+ 1) , fj=⇒ fj− aM, jV2(IΔr) , (I = 1, 2, ..., M − N; J = U1+ U2).

After the replacement we remove the relevant terms of a1, jxj+1, bM, jxj−M, and aM, jxj+Mfrom (7). Also we remove b1, jxj−1for I = M and J = U1+ 1, U1+ 2, ..., U1+ U2 owing to (2b). Now (7) becomes bL, jxj−L+ b1, jxj−1+ djxj+ a1, jxj+1+ aM, jxj+M= fj, (8a) ( j= N1+ 1, ..., N1+ L), b1, jxj₋₁+ djxj+ a1, jxj₊₁+ aM, jxj_+M= fj, (8b) ( j= N1+ L + 1, ..., N1+ M), bM, jxj−M+ b1, jxj−1+ djxj+ a1, jxj+1+ aM, jxj+M= fj, (8c) ( j= N1+ M + 1, ..., N1+ (U2− 1)M + N), bM, jxj−M+ b1, jxj−1+ djxj+ a1, jxj+1= fj, (8d) ( j= N1+ (U2− 1)M + N + 1, ..., N1+ N2).

Note that the first term of (8a) has this form because the grid point j− L is in the lower region. Regard dj, b1, j, and fjrelevant to the boundary conditions as the quantities after correction. Also we have left certain null a1, j’s and b1, j’s in (8).

In the upper region a general form is

bN, jxj−N+ b1, jxj−1+ djxj+ a1, jxj+1+ aN, jxj+N= fj. (9) Equation (2h) requires the replacement

After the replacement we remove a1, jxj+1from (9). Also, (2b) and (2d) require

b1, jxj₋₁=⇒ 0, (I = M; J = U1+ U2+ 1, ..., U),

aN, jxj_+N=⇒ 0, (I = M, M − 1, ..., M − N + 1; J = U). The resulting system of equations is

bM, jxj_−M+ b1, jxj₋₁+ djxj+ a1, jxj₊₁+ aN, jxj_+N = fj, (10a) ( j= N1+ N2+ 1, ..., N1+ N2+ N), bN, jxj−N+ b1, jxj−1+ djxj+ a1, jxj+1+ aN, jxj+N= fj, (10b) ( j= N1+ N2+ N + 1, ..., N1+ N2+ (U3− 1)N), bN, jxj−N+ b1, jxj−1+ djxj+ a1, jxj+1= fj, (10c) ( j= N1+ N2+ (U3− 1)N + 1, ..., n).

Regard fjas the quantity after correction of the boundary condition. Also, we left null a1, j’s and b1, j’s in (10) for convenience’ sake.

The total system of equations consists of (6), (8), and (10). The coefficients matrix has the form in Fig. 2.

Let us follow the procedure described in Section 3 of the previous paper.(4)_{We eliminate}

x1, x2, ..., xn−1in turn from the system of equations. As before, the index k means the set of equations that does not include xk. Let j = k + l, where l = 1, 2, ..., L in the lower region,

l= 1, 2, ..., M in the middle region, and l = 1, 2...., N in the upper region. Here we give only the equations for l= 1, which are required in obtaining the solution. Let j be k + 1. Starting from (6a), we obtained

dkjxj+ a1, jxj₊₁+ L_−1 i=L−k ak_{i, j}xj_+i+ aL, jxj_+L= fkj, (1 ≤ k ≤ L − 2), (11a) dk_jxj+ L_−1 i₌₁ ak_{i, j}xj+i+ aL, jxj+L= fjk, (L − 1 ≤ k ≤ N1− 1), (11b) dkjxj+ L_−1 i₌₁ aki, jxj+i+ aM, jxj+M= fkj, (k = N1), (11c)

L = 4

M = 11

N = 5

dkjxj+ s  i₌₁ aki, jxj+i+ M_−1 i=s aki, jxj+i+ aM, jxj+M = fjk, (N1+ 1 ≤ k ≤ N1+ L − 2), (11d) where s= N1+ L − k − 1 and s  = N1+ M − k. dkjxj+ a1, jxj₊₁+ M_−1 i_=s ak_{i, j}xj_+i+ aM, jxj_+M= fkj, (11e) (N1+ L − 1 ≤ k ≤ N1+ M − 2), dkjxj+ M_−1 i=1 ak_{i, j}xj_+i+ aM, jxj_+M= fjk, (11f) (N1+ M − 1 ≤ k ≤ N1+ N2+ N − M − 1), dkjxj+ s  i=1 aki, jxj_+i= fjk, (11g) (N1+ N2+ N − M ≤ k ≤ N1+ N2− 1), where s= N1+ N2+ N − k − 1. dkjxj+ N_−1 i=1 aki, jxj_+i+ aN, jxj_+N= fjk, (N1+ N2≤ k ≤ n − N − 1), (11h) dkjxj+ n_−k−1 i=1 aki, jxj_+i= fjk, (n − N ≤ k ≤ n − 2), (11i) dkjxj= fjk, (k = n − 1). (11j)

Let us address ourselves to dk j, a

k

i, j, bki, j, and fjk. For k= 1, we have

dkj= dj− c0al,k, (l = 1, 2, ..., L), (12a) aki, j= ai, j− c0al+i,k, (l = 1, 2, ..., L − 1; i = 1, 2, ..., L − l), (12b) bki, j= bi, j− c0al_−i,k, (l = 2, 3, ..., L; i = 1, 2, ..., l − 1), (12c) fkj = fj− c0fk, (l = 1, 2, ..., L), (12d) where c0_{= b} l, j/dk.

In case of k= 1 the required coefficients are dk jand f k j for l= 1 and L, a k L_{−1, j}for l= 1, and bk

L_{−1, j}for l= L. The required range of l changes in a complicated manner with increase of k. Obtaining the coefficients for all l’s greatly simplifies programming. The coefficients for k≥ 2 are given by dkj= d λ j− bμ_{l, j} dk−1 k aνl,k, (l = 1, 2, ..., m), (13a) aki, j= aλi, j− bμ_{l, j} dk−1 k aνl_+i,k, (l = 1, 2, ..., m − 1; i = 1, 2, ..., m − l) (13b) bk_{i, j}= bλ_{i, j}− b μ l, j d_kk−1a ν l−i,k, (l = 2, 3, ..., m; i = 1, 2, ..., l − 1), (13c) fk_j = f_jλ− b μ l, j dk−1 k f_kk−1, (l = 1, 2, ..., m), (13d)

where superscripts (λ, μ, ν) are k− 1 unless stated otherwise. Here replace m by L for 2 ≤ k ≤ N1, replace m by M for N1+ 1 ≤ k ≤ N1+ N2− 1, and replace m by N for N1+ N2≤ k ≤ n − 1. Exceptions appear on (λ, μ, ν) of (13).

 Exc.(l = 1): remove (λ, μ, ν) from equation of dk

j and remove (λ, μ) from equations of

ak

i, jand fjk.

 Exc.(i = l − 1): remove (λ, ν) from equation of bk i, j.

 Exc.(l = L), exc.(l = M), exc.(l = N): remove (λ, μ, ν) from equation of dk

j and remove (λ, μ) from equations of bk_{i, j}and f_jk.

 Exc.(i = L − l), exc.(i = M − l), exc.(i = N − l): remove (λ, ν) from equation of ak i, j. These exceptions are applied all together to the coefficients dk

j, a k

i, j, bki, j, and fjk. It is possible that two exceptions are applied simultaneously. Application of exceptions depends on the range of k. For 2 ≤ k ≤ L − 1, apply exc.(l = L), exc.(i = L − l), exc.(l = 1), and exc.(i = l − 1). For L ≤ k ≤ N1, apply exc.(l= L) and exc.(i = L − l). For N1+ 1 ≤ k ≤

N1+ L − 1, apply exc.(l = M) and exc.(i = M − l). For N1+ L ≤ k ≤ N1+ M − 1, apply exc.(l= M), exc.(i = M − l), exc.(l = 1), and exc.(i = l − 1). For N1+ M ≤ k ≤ N1+ N2, apply exc.(l= M) and exc.(i = M − l). For N1+ N2+ 1 ≤ k ≤ n − 1, apply exc.(l = N) and exc.(i= N − l). The procedure of obtaining the solution may be clear. Obtain the coefficients for k= 1, 2, ..., and we have xn, xn−1, ..., x1by using (11j) to (11a), and (6a).

Effectiveness of the method was examined by solving a sample problem. The choice of the reactor dimension is R = 200, R1 = 120, R2 = 100, H = 280, H1 = 100, and H2 = 200. The right -hand side of the system of (6), (8), and (10) was determined by substituting (x1, x2, ..., xn) of xj= cos  πIΔr 2R  sin  πJΔz H  (14)

into the left-hand side. The solution (X1, X2, ..., Xn) was first obtained forΔr = Δz = 0.5, for which n= 143041. The computer execution time t was 0.22s and the relative error ε defined by

ε =X − x x . was 3.8× 10−12.

The solution was also obtained forΔr = 0.5 and Δz = 0.2, for which n = 358201. The time t was 0.56s and the error ε was 1.1x10−11. The computer time and the error of this order are more than enough for PIC/MC simulations of processing plasmas.(1),(2)

The computer execution time for the present method was compared with that for SOR with Chebyshev acceleration.(6)This SOR is adopted because the optimal over relaxation pa-rameter is given explicitly as a function of Jacobi spectral radius. That is, there is no ambiguity in choosing the optimal relaxation parameter. The error limit ε of SOR is set to 10−7. In case of n= 143041 the computer time for SOR was 39 times larger and in case of n = 358201 the time for SOR is 92 times larger than that for the present method.

3. Three-dimensional plasma reactor

Three-dimensional PIC/MC simulations have been requested in plasma processing.(3),(7) A quick solver of 3D Poisson’s equation for electrical potentialΦ

∂2_Φ ∂x2 + ∂2_Φ ∂y2 + ∂2_Φ ∂z2 = g(x, y, z) (15)

plays a key role in plasma simulations. The domain is 0≤ x ≤ a, 0 ≤ y ≤ b, and 0 ≤ z ≤ c. Although we can choose arbitrary boundary conditions in the present method, we set for simplicityΦ = 0 on all boundary surfaces. Let us make a grid with size of Δx = a/(m + 1),Δy = b/(my+ 1), and Δz = c/(mz+ 1). If we think of ΦI,J,K asΦ (IΔx, JΔy, KΔz), the

central-difference representation of (15) becomes βx(ΦI_−1,J,K+ ΦI_+1,J,K)− 2  1+ βx+ βy  ΦI,J,K + βy(ΦI,J−1,K+ ΦI,J+1,K)+ (ΦI,J,K−1+ ΦI,J,K+1) = fI,J,K, (16)

where βx = (Δz/Δx)2, βy = (Δz/Δy)2, fI,J,K = (Δz)2g(IΔx, JΔy, KΔz), I = 1, 2, ..., m, J = 1, 2, ..., my, and K = 1, 2, ..., mz.

We introduce the serial grid number j instead of (I, J, K)

j= I + m(J − 1) + M(K − 1), (17)

where j= 1, 2, ..., n and M = mmy, Mmz(= n) being the number of unknowns. If we replace ΦI,J,K by xj, we have xj₊₁ = ΦI_+1,J,K, xj₋₁ = ΦI_−1,J,K, xj_+m = ΦI_,J+1,K, xj_−m = ΦI_,J−1,K,

xj+M = ΦI,J,K+1, and xj−M = ΦI,J,K−1. Then (16) takes the form

bM, jxj−M+ bm, jxj−m+ b1, jxj−1+ djxj+ a1, jxj+1+ am, jxj+m+ aM, jxj+M= fj, (18) where dj= −2(1 + βx+ βy), a1, j= b1, j = βx, am, j= bm, j= βy, aM, j= bM, j= 1.

Taking the null boundary conditions into consideration, we have, from (18) djxj+ a1, jxj₊₁+ am, jxj_+m+ aM, jxj_+M= fj, ( j = 1), (19a)

b1, jxj₋₁+ djxj+ a1, jxj₊₁+ am, jxj_+m+ aM, jxj_+M = fj, ( j = 2, 3, ..., m). (19b)

Here a1, jxj+1 = 0 for j = m from the boundary condition at x = a. But the term is left for convenience’ sake. For j > m we have

bm, jxj−m+ b1, jxj−1+ djxj+ a1, jxj+1+ am, jxj+m+ aM, jxj+M= fj, (19c) ( j= m + 1, m + 2, ..., M − m),

where b1, jxj₋₁ = 0 for j with I = 1, J = 2, 3, ..., my− 1, and K = 1. Also, a1, jxj₊₁ = 0 for j with I= m, J = 2, 3, ..., my− 1,and K = 1. For j > M − m we have

bm, jxj_−m+ b1, jxj₋₁+ djxj+ a1, jxj₊₁+ aM, jxj_+M = fj, (19d) ( j= M − m + 1, ..., M),

where b1, jxj−1= 0 for j = M − m + 1 and a1, jxj+1= 0 for j = M. For j > M we have

bM, jxj_−M+ bm, jxj_−m+ b1, jxj₋₁+ djxj

+ a1, jxj+1+ am, jxj+m+ aM, jxj+M = fj, (19e) ( j= M + 1, M + 2, ..., n − M),

where b1, jxj₋₁ = 0 for j with I = 1, a1, jxj₊₁ = 0 for j with I = m, bm, jxj_−m = 0 for j with

J= 1, and am, jxj_+m= 0 for j with J = my. For j > n − M + 1 we have

bM, jxj−M+ bm, jxj−m+ b1, jxj−1+ djxj+ a1, jxj+1+ am, jxj+m= fj, (19f) ( j= n − M + 1, n − M + 2, ..., n),

m = 5

M = 25

Fig. 3 Coefficients matrix for 3D reactor

where b1, jxj−1 = 0 for j with I = 1, a1, jxj+1 = 0 for j with I = m, bm, jxj−m = 0 for j with

J= 1, and am, jxj+m= 0 for j with J = my.The coefficients matrix of (19) is shown in Fig.3. It has nonzero elements on seven diagonal lines.

Let us proceed as usual. Obtain x1from (19a) and substitute it into equations for j = 2, m+ 1, and M + 1, and we have

dkjx2+ a1, jxj₊₁+ akm_{−1, j}xj_+m−1+ am, jxj_+m + ak M−1, jxj+M−1+ aM, jxj+M = fjk, (l = 1), (20a) bkm−1, jx2+ b1, jxj−1+ dkjxj+ a1, jxj+1+ am, jxj+m + ak M_{−m, j}xj+M−m+ aM, jxj+M = fkj, (l = m), (20b) bk_{M−1, j}x2+ bkM−m, jxj−M+m+ bm, jxj−m+ b1, jxj−1+ dkjxj + a1, jxj₊₁+ am, jxj_+m+ aM, jxj_+M = fkj, (l = M), (20c) where j = k + l and k = 1.The coefficients dk

j, a k

i, j, bki, j, and fjkfor k= 1 are given by (12), where L is replaced by M. The set of k= 1 includes also (19b) for j = 3, (19c) for j = m + 2, and (19e) for j= M + 2. If we eliminate x2from the set of k= 1, we have the set of k = 2. Here we give only the equations of l= 1 for all sets, j being k + 1.

dkjxj+ a1, jxj₊₁+ m_−1 i=m−k ak_{i, j}xj_+i+ am, jxj_+m+ M_−1 i=M−k ak_{i, j}xj_+i + aM, jxj+M = fjk, (1 ≤ k ≤ m − 2), (21a) dkjxj+ m_−1 i=1 aki, jxj_+i+ am, jxj_+m+ M_−1 i=M−k aki, jxj_+i + aM, jxj+M= fjk, (m − 1 ≤ k ≤ M − m − 1), (21b) dkjxj+ M_−1 i=1 ak_{i, j}xj+i+ aM, jxj+M= fjk, (M − m ≤ k ≤ n − M − 1), (21c)

dkjxj+ n−k−1_

i₌₁

aki, jxj+i= fjk, (n − M ≤ k ≤ n − 2), (21d)

dkjxj= fjk, (k = n − 1). (21e)

The coefficients for k ≥ 2 are given by (13), where m is to be replaced by M. Exceptions on the removal of (λ, μ, ν) are as follows. Exc.(l= 1) is the same as in Section 2. Other exceptions are summarized as follows.

 Exc.(l = m): remove (λ, μ, ν) from equation of dk

j and remove (λ, μ) from equations of

ak

i, j, bki, jand fkj.

 Exc.(l = M): remove (λ, μ, ν) from equation of dk

j and remove (λ, μ) from equations of bk

i, jand fjk.

 Exc.(i = m − l), exc.(i = M − l): remove (λ, ν) from equation of ak i, j.  Exc.(i = l − 1), exc.(i = l − m): remove (λ, ν) from equation of bk

i, j.

As before, in certain cases two exceptions are applied simultaneously. Application of exceptions depends on the range of k . Apply all exceptions for 2 ≤ k ≤ m − 1, apply exc.(l= m), exc.(l = M), exc.(i = m − l), exc.(i = M−l), and exc.(i = l−m) for m ≤ k ≤ M−m, apply exc.(l= M) and exc.(i = M − l) for M − m + 1 ≤ k ≤ n − M, and apply no exceptions for n− M + 1 ≤ k ≤ n − 1. The coefficients dk_j, ak_{i, j}, bk_{i, j}and f_jkfor l= 1, 2, ..., M are required for k≤ n − M − 1 and those for l = 1, 2, ..., n − k are required for k ≥ n − M.

Validity and effectiveness of the present method were examined by a sample calculation. Substituting xj= sin  πIΔx a  sin  πJΔy b  sin  πKΔz c  (22)

into the left-hand side of the system of (19), we determined the right-hand side. Let (X1, X2, ..., Xn) be the obtained solution. The relative error is defined by as before. We also obtained the solution by using the FFT method; (16) are doubly Fourier-sine transformed with respect to I and J, and then the resulting difference equation is solved by use of the Thomas algorithm. The ASL library(8)was used in applying FFT. The reason of not using SOR this time is in that we could not find 3D SOR scheme whose optimal over relaxation parameter is given explicitly. The computation was carried out for m = my = mz = 50, i.e. n = 125000. The computer execution time t and the relative error ε for the two methods were compared. The time t was nearly the same for the two methods and the error ε for the FFT was 10−14, one order smaller than that for the present method. Note ,however, that the FFT method is applicable only when the boundary conditions are very simple.

4. Conclusion

The direct method developed in the first paper(4) _{to solve the finite-difference systems} for Poisson’s equation is extended to the cases with more complex boundary conditions. The first case is the axisymmetrical plasma reactor (GEC cell) of the international standard and the second is the three dimensional plasma reactor. In the first case the computer execution time for the present method is 1 to 3% of the time for SOR with Chebyshev acceleration. In the second case the finite-difference system is solved by the present method and also by the FFT method. The computer time for the present method is nearly the same as that for FFT method. Note, however, that the FFT method is applicable only when boundary conditions are very simple.

Acknowledgements

The author would like to express his thanks to Prof. S. Yamamoto in Tohoku University for accepting him as a guest scientist. He also expresses his thanks to Dr. Y. Sasao for advice on programming. All computations were performed by use of the vector computer NEC SX-9 in Cyber Science Center, Tohoku University. The author thanks to Mr. T. Yamashita for his kind advice on a speedup technique of SX-9.

References

( 1 ) Birdsall, C. K. and Langdon, A. B., Plasma Physics via Computer Simulation, McGraw-Hill, New York, 1985.

( 2 ) Nanbu, K., Probability theory of electron-molecule, ion-molecule, molecule-molecule, and Coulomb collisions for particle modeling of materials processing plasmas and gases, IEEE Trans. Plasma Sci. 28 (2000) 971-990.

( 3 ) Kadlec, S., Computer simulation of magnetron sputtering - Experience from the indus-try, Surf. Coat. Tech. 202 (2007) 895-903.

( 4 ) Nanbu, K., Direct solutions of finite-difference systems for Poisson’s equation I. Simple cases, J. Comp. Sci. Tech. 6 (2012) 147-156.

( 5 ) Lymberopoulos, D. P. and Economou, D. J., Two-dimensional self-consistent radio fre-quency plasma simulations relevant to the gaseous electronics conference rf reference cell, J. Res. Nat. Inst. Stand. Technol. 100 (1995) 473-494.

( 6 ) Numerical Recipes in Fortran 77: The Art of Scientific Computing, Cambridge Univer-sity Press, (1992), 857-860.

( 7 ) Nanbu, K. and Ohshita, T., Synthetic simulation of plasma formation, target erosion, and film deposition in a large magnetron sputtering apparatus, Vacuum 87 (2013) 103-108. ( 8 ) Advanced Scientific Library (ASL/SX), NEC Corporation, 7th ed., 2002.