A conservative parallel difference method for 2-dimension diffusion equation

Accepted Manuscript A conservative parallel difference method for 2-dimension diffusion equation Dongxu Jia, Zhiqiang Sheng, Guangwei Yuan

PII: DOI: Reference:

S0893-9659(17)30330-0 https://doi.org/10.1016/j.aml.2017.11.004 AML 5366

To appear in:

Applied Mathematics Letters

Received date : 2 November 2017 Accepted date : 8 November 2017 Please cite this article as: D. Jia, Z. Sheng, G. Yuan, A conservative parallel difference method for 2-dimension diffusion equation, Appl. Math. Lett. (2017), https://doi.org/10.1016/j.aml.2017.11.004 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Manuscript Click here to view linked References

A conservative parallel difference method for 2-dimension diffusion equation Dongxu Jiaa,, Zhiqiang Shengb , Guangwei Yuanb,∗ b Laboratory

a The Graduate School of China Academy of Engineering Physics, P.O.Box 2101, Beijing, China, 100088 of Computational Physics, Institute of Applied Physics and Computational Mathematics, P.O. Box 8009, Beijing, China, 100088

Abstract In this paper, a conservative parallel difference scheme, which is based on domain decomposition method, for 2-dimension diffusion equation is proposed. In the construction of this scheme, we use the numerical solution on the previous time step to give a weighted approximation of the numerical flux. Then the sub-problems with Neumann boundary are computed by fully implicit scheme. What’s more, only local message communication is needed in the program. We use the method of discrete functional analysis to give the proof of the unconditional stability and second-order convergence accuracy. Some numerical tests are given to verify the theory results. Keywords: parallel difference, diffusion equation, conservative, domain decomposition

1. Introduction Time-dependent diffusion equations are widely used in many fields of science and engineering. To solve diffusion equations efficiently, the classic fully implicit method is usually a natural choice because of its unconditional stability [1]. However a global linear algebraic systems resulting from the implicit method must be solved at each time level which cannot be implemented directly and naturally on parallel computers. On the other hand, the classic explicit scheme can not be used in some applications due to its severe time step-length restriction though it has intrinsic parallelism. Note that conservation is also an important property for diffusion equation, so it is reasonable to devise a conservative scheme for it. In summary it is necessary to devise new schemes for diffusion equations satisfying the following properties: (i) unconditional stability; (ii) second-order accuracy in space; (iii) intrinsic parallelism; (iv) conservation. Several early attempts are carried out in [2–8]. Especially in [9, 10] some conservative schemes are devised and numerical tests show that they have good stability and second-order accuracy, but no theoretical proof is given. In [11] a conservative parallel scheme for one-dimensional diffusion equation is devised and proved theoretically to be of unconditional stability and second-order accuracy. But, as we know, for any conservation parallel scheme of multidimensional diffusion problems there is few theoretical result up to now. For two-dimensional problems, another approach based on splitting technique is proposed in [12] to devise parallelization procedures, which are different from all the parallel difference schemes mentioned above. In this work, we construct a new parallel scheme for diffusion equations by predicting a numerical flux on the interface with the value on previous time-level. Then the fully implicit scheme is carried out on the sub-domains. This method can easily be extended to high-dimensional problems ( for example 2 dimension ones). Both numerical tests and theoretical analysis demonstrate that the scheme is conservative, unconditionally stable and has second-order accuracy. This article is organized as follows. In section 2, a conservative parallel difference scheme is given for 2-dimension and we also list main results in this part. In section 3, both unconditional stability and convergence are demonstrated for 1-dimension problem. We extend the proof to 2-dimension in section 4. In section 5 some numerical tests are given to verify the theoretical results. We conclude the paper in section 6. ∗ Corresponding

author Email address: [email protected] (Guangwei Yuan)

Preprint submitted to Elsevier

November 1, 2017

2. Some Notations, Parallel Schemes and Main Results In this paper, we consider the following model problem Ut (x, t) + ∇ · Q(x, t) = 0, (x, t) ∈ Ω × (0, T ],

(2.1)

Q(x, t) · n = 0, (x, t) ∈ ∂Ω × (0, T ],

(2.3)

∇U(x, t) + Q(x, t) = 0, (x, t) ∈ Ω × (0, T ],

(2.2)

U(x, 0) = U0 (x), x ∈ Ω,

(2.4)

where Ω ⊂ R2 is a bounded domain and n is the outward normal direction on ∂Ω. For simplicity, we take Ω = [0, 1]2 throughout this work. Let (x1 , x2 ) be a point in R2 . We divide the domain Ω with J1 segments in x1 direction and x2 direction with J2 ones, then we set h1 = J11 , h2 = J12 be the space step-lengths. We generate the mesh on Ω by the mesh points (x j1 , x j2 ) = ( j1 h1 , j2 h2 ), where j1 = 0, 1, 2, · · · , J1 ; j2 = 0, 1, 2, · · · , J2 . For a function F(x1 , x2 ), let F j1 , j2 = F(x j1 , x j2 ) be the mesh function, then we can define some difference operators as follows ∆1+ F j1 , j2 = F j1 +1, j2 − F j1 , j2 , ∆2+ F j1 , j2 = F j1 , j2 +1 − F j1 , j2 , where jk = 21 , 1, 32 , 2, 52 , · · · , Jk − 12 , (k = 1, 2). And ∆1− F j1 , j2 = F j1 , j2 − F j1 −1, j2 , ∆2− F j1 , j2 = F j1 , j2 − F j1 , j2 −1 , where jk = 23 , 2, 52 , · · · , Jk − 21 , Jk , Jk + 12 , (k = 1, 2). Let τ = T/Nt , tn = nτ, where Nt is a positive integer. Then for a function F(t), denote F n = F(tn ), (n = 0, 1, 2, · · · , Nt ). Let λ1 = τ/h21 , λ2 = τ/h22 be the mesh ratios. Now we give the conservative parallel difference scheme for 2-dimension problem. For the sake of simplicity, we separate the domain with a cross and assume that the separating interface on x1 , x2 directions are j∗1 h1 , j∗2 h2 , respectively. Where we suppose that j∗1 , j∗2 are integers with 3 < j∗1 < J1 − 3 and 3 < j∗2 < J2 − 3 . Then the 2-dimension scheme can be written as follows n un+1 j1 , j2 − u j1 , j2

∆1− qn+1 j +1,j 1

2

2

! ∆2− qn+1 j ,j +1 1

2

2

, ( j1 = 1, 2, ..., J1 , j2 = 1, 2, ..., J2 ),

(2.5)

=

−

qn+1 j +1,j

=

−

qn+1 j ,j +1

=

−

qn+1 j∗ + 1 , j

=

ω1

qn+1 j , j∗ + 1

=

ω2

qn+1 1 ,j

=

qn+1 j ,1

qn+1 = 0, ( j2 = 1, 2, · · · , J2 ), J +1,j

=

(2.11)

u0j1 , j2

qn+1 = 0, ( j1 = 1, 2, · · · , J1 ), j ,J + 1

=

U0 ( j1 h1 , j2 , h2 ), ( j1 = 1, 2, ..., J1 ; j2 = 1, 2, ..., J2 ),

(2.12)

τ

1

1

1

1

2

2

2

2

2

2

2

2

2

2

1 2

h1 ∆1+ un+1 j1 , j2 h1 ∆2+ un+1 j1 , j2

1

1

+

h2

, ( j1 , j∗1 , j2 = 1, 2, ..., J2 ),

, ( j2 , j∗2 , j1 = 1, 2, ..., J1 ), h2 qnj∗ + 3 , j + qnj∗ − 1 , j ! ∆1+ unj1 , j2 ! 2 1 2 2 − + (1 − ω1 ) 1 2 , ( j2 = 1, 2, ..., J2 ), h1 2 qnj , j∗ + 3 + qnj , j∗ − 1 ! ∆2+ unj1 , j2 ! 1 2 2 1 2 2 − + (1 − ω2 ) , ( j1 = 1, 2, ..., J1 ), h2 2 2

2

2

2

where ωk are weights depending on the mesh ratio λk = τ/h2k , (k = 1, 2).

2

(2.6) (2.7) (2.8) (2.9) (2.10)

Remark 2.1. To devise a conservative scheme we should define all the unknowns in the center of each mesh element. That is, we put all the unknowns unj1 , j2 (n = 0, 1, 2, · · · ) on the following locations,   ( j1 − 1/2)h1 , ( j2 − 1/2)h2 ,

(2.13)

where j1 = 1, 2, · · · , J1 ; j2 = 1, 2, · · · , J2 . Then it is easy to see that all the numerical fluxes qn (n = 0, 1, 2, · · · ) are defined on the faces of each mesh element, and the flux locations should also be calculated by (2.13). Theorem 1. The scheme (2.5)-(2.12) is unconditionally stable with L2 discrete norm [13] when ωk satisfy (   1/(2λk + 5), 1 , when λk < 1   , k = 1, 2. ωk ∈ 1/(2λk + 5), 1/(2λk − 1) , when λk ≥ 1

Theorem 2. The scheme (2.5)-(2.12) has second-order accuracy with L2 discrete norm when ωk satisfy (   1/(2λk + 5 − 4), 1 , when λk < 1 −    , k = 1, 2, ωk ∈ 1/(2λk + 5 − 4), 1/(2λk − 1 + 2) , when λk ≥ 1 −  where  ∈ (0, 1] is a small positive constant.

Corollary 1. Scheme (2.5)-(2.12) is unconditionally stable and has second-order accuracy when ωk = 1, 2) and this choice coincides exactly with the one dimension method in [11].

1 1+2λk ,

(k =

By examining the scheme for 2-dimension problem, it is obvious that the following theorem holds. P P P P Theorem 3. Scheme (2.5)- (2.12) is a conservative scheme, that is Jj11=1 Jj22=1 unj1 , j2 h1 h2 = Jj11=1 Jj22=1 u0j1 , j2 h1 h2 holds, where n = 1, 2, · · · , Nt . 3. Theoretical Analysis in One Dimension 3.1. Proof of unconditional stability In this section, we give the proof of stability for the scheme. And we will lose nothing interest by confining our attention to one dimension case. We omit all the subscripts having relationship with the dimension, then we denote h for the space step-length and so on. Rewrite the scheme in one dimension with two sub-domains as follows (we take the separating interface as kh): τ−1 (un+1 − unj ) = −h−1 ∆− qn+1 , ( j = 1, 2, · · · , J), j j+ 1 2

qn+1 = −h−1 ∆+ un+1 j , ( j = 1, 2, · · · , k − 1, k + 1, · · · , J), j+ 1 2

n n  qk+ 32 + qk− 21   ∆+ unk  + (1 − ω) , qn+1 = ω − k+ 12 h 2 qn+1 = qn+1 = 0, 1 J+ 1 2

(3.1) (3.2) (3.3) (3.4)

2

u0j = U0 (x j ), j = 1, 2 · · · , J.

(3.5)

where the weight ω is to be determined in the stability proof. n+1 Before we give the proof of the stability, we rewrite the numerical flux qk+ 1 to 2

qn+1 k+ 12

n n n n  ∆+ unk   qk+ 32 + qk− 12   ∆+ unk   qk+ 23 + qk− 12  ∆+ un+1 ∆+ un+1 k k =ω − + (1 − ω) = −ω +ω +ω − + (1 − ω) h 2 h h h 2 n+1     ∆+ u k 1−ω n = −ω − ωλ∆+ qn+1 − qn+1 + qk+ 3 + qnk− 1 . (3.6) k+ 12 k− 21 2 2 h 2

3

Proof. Multiply at both side of the scheme (3.1) by un+1 j h and sum over with j = 1, 2, · · · , J, to get J

J

1X 1 X n+1 (u j − unj )un+1 ∆− qn+1 un+1 h j h+ j+ 12 j τ j=1 h j=1

0=

J

J

J

J

X X i X 1 h X n+1 2 (unj )2 h + (un+1 − unj )2 h + ∆− qn+1 (u j ) h − un+1 . j j+ 12 j 2τ j=1 j=1 j=1 j=1

=

(3.7)

Now we deal with the last item on the right-hand side in (3.7). When we notice (3.6), then one gets J X j=1

∆− qn+1 un+1 = j+ 1 j 2

=

J−1 X

j=1, j,k J−1 X j=1

(qn+1 )2 h − qn+1 ∆ un+1 j+ 1 k+ 1 + k 2

2

(qn+1 )2 h + j+ 1 2

   1 − ω n+1 2 1 − ω n n+1 n+1 (qk+ 1 ) h + λ∆+ qn+1 qk+ 3 + qnk− 1 qn+1 h. 1 − q 1 q 1h − k+ k− k+ k+ 12 2 2 2 2 2 2 ω 2ω

Using the Cauchy inequality 2ab ≤ a2 + b2 , we obtain J X j=1

un+1 ≥ ∆− qn+1 j+ 1 j 2

J−1 X j=1

− Then with the identity (a + b)c =

)2 h + (qn+1 j+ 1 2

  1 − ω n+1 2 n+1 n+1 (qk+ 1 ) h + λh qn+1 − 2λh(qn+1 )2 3 + q 1 q k+ k− k+ 21 k+ 12 2 2 2 ω

 (1 − ω)h  n 2 2 (qk+ 3 ) + (qnk− 1 )2 + 2(qn+1 . 1) k+ 2 2 2 4ω

a2 +b2 +2c2 2

−

(a−c)2 +(b−c)2 , 2

and (3.7), we can obtain

J J   i 1 X n+1 2 1 X n2 (1 − ω) h n+1 2 2 (u j ) h − (u j ) h + (qk+ 3 ) + (qn+1 − (qnk+ 3 )2 + (qnk− 1 )2 h 1) k− 2 2 2 2 2τ j=1 2τ j=1 4ω X h1 − ω i h λ (1 − ω) i  2 2 n+1 2 )2 h − − λ + 1 (qn+1 − + 1 (qn+1 ≤− (qn+1 h 1) h − 3 ) + (q 1) k+ k+ k− j+ 21 2 2 2 2ω 2 4ω j,k,k±1

−

J  1 X n+1 λh  n+1 2 n+1 n+1 2 (u j − unj )2 h + (qk+ 3 − qn+1 . 1 ) + (q 1 − q 1) k+ k− k+ 2 2 2 2 2τ j=1 2

(3.8)

By using the parallel scheme again, the sign of the last three items in (3.8) is nonpositive. So if the following relationship holds: 1−ω λ (1 − ω) − λ + 1 ≥ 0; − + 1 ≥ 0; 0 ≤ ω ≤ 1, 2ω 2 4ω which mean that ω satisfies the following conditions: (   1/(2λ + 5), 1 , when λ < 1  ω∈  , 1/(2λ + 5), 1/(2λ − 1) , when λ ≥ 1

we can obtain the following inequality from (3.8), J

J

  i 1 X n+1 2 1 X n2 (1 − ω) h n+1 2 2 n 2 n 2 (u j ) h − (u j ) h + (qk+ 3 ) + (qn+1 − (q h ≤ 0. 1) 3 ) + (q 1) k− 2 k+ 2 k− 2 2 2τ j=1 2τ j=1 4ω

Hence, the parallel difference scheme is unconditionally stable under the following energy norm J

En+1 :=

 1 X n+1 2 (1 − ω)τ  n+1 2 2 (u j ) h + (qk+ 3 ) h + (qn+1 1) h , k− 2 2 2 j=1 4ω

and the norm is an equivalent norm with L2 discrete norm. This completes the proof. 4

3.2. Proof of second-order convergence Let qnj+ 1 be the numerical flux and the exact flux is Qnj+ 1 . Denote ηnj+ 1 = Qnj+ 1 − qnj+ 1 . Then the error equations 2

2

can be written as

2

2

2

τ−1 (ξn+1 − ξnj ) = −h−1 ∆− ηn+1 + T n+1 j j , j = 1, 2, · · · , J, j+ 1

(3.9)

2

ηn+1 = −h−1 ∆+ ξn+1 + S n+1 , j , k, j j+ 1 j+ 1 2

2

ηn+1 =ω − k+ 1 2

! ∆+ ξkn + (1 − ω) h

ηnk+ 3 2

(3.10) +

ηnk− 1 2

2

= ηn+1 = 0, ξ0j = 0. ηn+1 1 J+ 1 2

2

!

n+1 + Sg , k+ 1

(3.11)

2

(3.12)

to We can rewrite the expression of ηn+1 k+ 1 2

! ηnk+ 3 ∆+ ξkn+1 2 − ωλ ∆+ ∆− ηn+1 + (1 − ω) k+ 12 2 h ! ηnk+ 3 ∆+ ξkn+1 2 n+1 =: − ω − ωλ ∆+ ∆− ηk+ 1 + (1 − ω) 2 h

ηn+1 =−ω k+ 1

+ ηnk− 1 ! 2

2 + ηnk− 1 ! 2

2

τ n+1 + Sg + ω ∆+ T kn+1 k+ 12 h n+1 + ωS k+ 1, 2

= O(τ + h2 ), T n+1 = O(τ + h2 ), j = 1, 2, · · · , J − 1. and it is obvious that S n+1 j j+ 1 2

With the similar method used in stability proof, we can get the convergence result.

4. Theoretical Analysis in 2-Dimension Multiply at both side of the scheme (2.5) by un+1 j1 , j2 h1 h2 , and sum over for each dimension, to obtain ! J2 J1 X J2 J1 X X 1 1X 1 n n+1 n+1 n+1 n+1 (un+1 − u )u h h = − + ∆ q ∆ q 1− j + 1 , j 2− j , j + 1 u j1 , j2 h1 h2 j1 , j2 j1 , j2 1 2 1 2 2 1 2 2 τ j =1 j =1 j1 , j2 h h 1 2 j1 =1 j2 =1 2 1 ! # ! # " " J J J J 2 2 1 1 X X X X 1 1 n+1 n+1 n+1 u h h + u h = − − ∆1− qn+1 ∆ q 2− j , j + 1 j1 , j2 1 2 j1 , j2 2 h1 . j1 + 12 , j2 1 2 2 h1 h2 j =1 j =1 j =1 j =1 2

1

1

2

Using the result in one dimension, we obtain En+1 ≤ En , (n = 0, 1, 2, · · · , Nt ), where

J1 X J2 1X (un+1 )2 h1 h2 2 j =1 j =1 j1 , j2 1 2 J1 " J2 " X X # # (1 − ω2 )τ  n+1 (1 − ω1 )τ  n+1 2 n+1 2 2 n+1 2 + (q j , j∗ + 3 ) + (q j , j∗ − 1 ) h1 h2 + (q j∗ + 3 , j ) + (q j∗ − 1 , j ) h1 h2 , 1 2 2 1 2 2 1 2 2 1 2 2 4ω2 4ω1 j =1 j =1

En+1 :=

1

2

and ωk satisfies the conditions listed in theorem 1, So we complete the proof in 2-dimension. The proof of convergence in 2-dimension is similar to the one dimension case . 5. Numerical Experiment

In this part, some numerical tests are given to verify the theoretical results in this work. We take the following two dimension model problem to test the 2D algorithm Ut − U xx − Uyy = (2π2 − A) exp(−At) sin πx sin πy, (x, y, t) ∈ (0, 1)2 × (0, T ], U x (0, y, t) = −U x (1, y, t) = π exp(−At) sin πy, (y, t) ∈ (0, 1) × (0, T ],

Uy (x, 0, t) = −Uy (x, 1, t) = π exp(−At) sin πx, (x, t) ∈ (0, 1) × (0, T ],

U(x, y, 0) = sin πx sin πy, (x, y) ∈ (0, 1)2 , 5

where A = 0.001 and the model problem has the exact solution U(x, y, t) = exp(−At) sin πx sin πy. We divide the domain Ω = [0, 1]2 into 4 × 4 sub-domains with two sets of orthogonal equidistant parallel lines , which parallel to x-axis or y-axis. At each sub-domain we choose N = N x = Ny = 8, 16, 32, 64, and let h = h x = hy = 1 to test the algorithm in two dimension. 1/N, λ = λ x = λy = τ/h2 and T = Nt τ. In the scheme we take ω1 = ω2 = 1+2λ All the results are demonstrated in the table 5.1 and table 5.2. It is obvious that the scheme is unconditionally stable and has second-order accuracy. Table 5.1: Accuracy for 2d: λ = 1, T = 1.

N ku − UkL2 Ratio

8 3.2482E-03

λ ku − UkL2

1 1.8618E-05

16 8.0689E-04 4.0256

32 2.0128E-04 4.0087

Table 5.2: Stability for 2d: N = 32, Nt = 1000.

10 1.3609E-04

100 1.3152E-03

64 5.0279E-05 4.0033

1000 1.2412E-02

6. Conclusion Conservation is an important property for parabolic equations. In this work, we propose a conservative parallel difference scheme for 2-dimension diffusion equation. The theoretical analysis is also given by using the skills of discrete functional analysis. Some numerical tests are presented to verify the theoretical results. It is our ongoing work to extend the scheme on distorted meshes. 7. Acknowledgements This work is partially supported by the National Natural Science Foundation of China (11571047, 11571048), NASF (U1630249), Science Challenge Project ( No.JCKY2016212A502 ) and the Foundation of CAEP (2015B0202042). References [1] K. W. Morton, D. F. Mayers, Numerical Solution of Partial Differential Equations: An Introduction, Cambridge University Press, 2005. [2] C. N. Dawson, Q. Du, T. F. Dupont, A finite difference domain decomposition algorithm for numerical solution of the heat equation, Mathematics of Computation 57 (195) (1991) 63–71. [3] C. N. Dawson, T. F. Dupont, Explicit/implicit, conservative domain decomposition procedures for parabolic problems based on block-centered finite differences, SIAM Journal on Numerical Analysis 31 (4) (1994) 1045–1061. [4] G. Yuan, F. Zuo, Parallel difference schemes for heat conduction equations., International Journal of Computer Mathematics 80 (8) (2003) 993–997. [5] G. Yuan, L. Shen, Stability and convergence of the explicit-implicit conservative domain decomposition procedure for parabolic problems, Computers & Mathematics with Applications 47 (4) (2004) 793–801. [6] H. L. Liao, H. S. Shi, Z. Z. Sun, Corrected explicit-implicit domain decomposition algorithms for two-dimensional semilinear parabolic equations, Science China Mathematics 52 (11) (2009) 2362. [7] H. Shi, H. Liao, Unconditional stability of corrected explicit/implicit domain decomposition algorithms for parallel approximation of heat equations, SIAM Journal on Numerical Analysis 44 (4) (2006) 1584–1611. [8] Z. Sheng, G. Yuan, X. Hang, Unconditional stability of parallel difference schemes with second order accuracy for parabolic equation, Applied Mathematics & Computation 184 (2) (2007) 1015–1031. [9] G. Yuan, X. Hang, Conservative parallel schemes for diffusion equations, Jisuan Wuli/Chinese Journal of Computational Physics 27 (4) (2010) 475–491. [10] G. Yuan, Y. Yao, L. Yin, A conservative domain decomposition procedure for nonlinear diffusion problems on arbitrary quadrilateral grids, SIAM Journal on Scientific Computing 33 (3) (2011) 1352–1368. [11] S. Zhu, Conservative domain decomposition procedure with unconditional stability and second-order accuracy, Applied Mathematics & Computation 216 (11) (2010) 3275–3282. [12] Z. Zhou, D. Liang, The mass-preserving and modified-upwind splitting DDM scheme for time-dependent convection-diffusion equations, Journal of Computational and Applied Mathematics 317 (2017) 247–273. [13] Y. Zhou, Applications of Discrete Functional Analysis to the Finite Difference Method, International Academic Publishers, 1990.

6