# A New Efficient Algorithm for the 2D WLP-FDTD Method Based on Domain Decomposition Technique

Authors: Bo-Ao Xu; Yan-Tao Duan; Bin Chen; Yun Yi; Kang Luo
Journal: International Journal of Antennas and Propagation
 (2016)

---

## Abstract

This letter introduces a new efficient algorithm for the two-dimensional weighted Laguerre polynomials finite difference time-domain (WLP-FDTD) method based on domain decomposition scheme. By using the domain decomposition finite difference technique, the whole computational domain is decomposed into several subdomains. The conventional WLP-FDTD and the efficient WLP-FDTD methods are, respectively, used to eliminate the splitting error and speed up the calculation in different subdomains. A joint calculation scheme is presented to reduce the amount of calculation. Through our work, the iteration is not essential to obtain the accurate results. Numerical example indicates that the efficiency and accuracy are improved compared with the efficient WLP-FDTD method.


---

## Body

## 1. Introduction

The Laguerre-based finite difference time-domain (WLP-FDTD) method [1] is widely used in solving electromagnetic problems with fine structures. As it is not constrained by the Courant-Friedrich-Levy condition, the WLP-FDTD method is unconditionally stable.However, the WLP-FDTD method often produces a huge sparse matrix equation in the process of solving the fine structure problem. Directly solving the associate matrix is a challenging work. In order to settle this problem, He et al. introduced the domain decomposition finite difference technique to the WLP-FDTD method [2]. The huge sparse matrix is decomposed into several small ones. Duan et al. investigated an efficient WLP-FDTD method [3], in which the huge sparse matrix is turned into tridiagonal matrixes by introducing a perturbation term. CPU time and memory storage can be saved greatly. However, the splitting error caused by the perturbation term turns to be pronounced in the fine structure area. To improve this efficient method, Chen et al. studied a new splitting scheme by introducing a new perturbation term [4]. However, the iteration is required in this method, which increases the computational time.In this letter, we introduce a new efficient algorithm for the 2D WLP-FDTD method. The whole computational domain is decomposed into several subdomains utilizing the domain decomposition scheme. In order to eliminate the splitting error, the conventional WLP-FDTD method [1] is applied to the subdomains where the fields have larger spatial derivatives. And to speed up the calculation, the efficient WLP-FDTD method [4] is used to the other subdomains where the splitting error is limited. Meanwhile, we devise a joint calculation scheme to further reduce the computational time. Through using this approach, the iteration is no longer necessary to get the accurate results. Numerical example shows that the accuracy can be improved and the CPU time is reduced to about 52.7% of the original one [4].


## 2. Theories

### 2.1. Domain Decomposition Technique

The domain decomposition technique is similar to that in [2]; the computational domain is decomposed into nine subdomains, as shown in Figure 1. However, in this paper, the conventional WLP-FDTD method is applied to D5 (fine structure exists). Meanwhile, the efficient WLP-FDTD method is used to the other subdomains and the interface Γ, where Γ=Γ1∪Γ2∪Γ3∪Γ4.Figure 1
Nine subdomains in the computational domain.In the simple and lossless medium, Maxwell’s equations of 2DTEz wave for the Laguerre-FDTD can be deduced as [4](1a)-abΔx2Eyqi-1,j+1+2abΔx2Eyqi,j-abΔx2Eyqi+1,j=2abΔxΔy∑k=0q-1Exki,j+1-Exki-1,j+1-Exki,j+Exki-1,j-2∑k=0q-1Eyki,j+2aΔx∑k=0q-1Hzki,j-Hzki-1,j-aJyqi,j(1b)-abΔy2Exqi,j-1+1+2abΔy2Exqi,j-abΔy2Exqi,j+1=-abΔxΔyEyqi+1,j-Eyqi+1,j-1-Eyqi,j+Eyqi,j-1-2∑k=0q-1Exki,j-2aΔy∑k=0q-1Hzki,j-Hzki,j-1(1c)a=2sε0,b=2sμ0,where ε0 is the free space electric permittivity, μ0 is the free space magnetic permeability, and Jyqi,j is the electric excitation source. The time-scale factor and order of the Laguerre polynomials are defined as s and q, respectively.Rewriting (1a) and (1b) as the matrix equations, we have(2a)Ayq=f1q-1(2b)Bxq=f2q-1,where vector x and y represent the unknown field variables Exqi,j and Eyqi,j, respectively.The conventional WLP-FDTD equations [1] can be written as(3)Cxyq=f3q-1,where vector xy represents the unknown field variables Eyqi,j and Exqi,j.According to (1a), in the subdomains D1~D3 and D7~D9, A in (2a) correspond to the interior to interior coupling. So firstly, yq in these subdomains can be solved by (2a) independently. Therefore one can have the sparse matrix system as(4)A440000A40A66000A600A2200A2000A880A80000A55A5A4∗A6∗A2∗A8∗A5∗AΓy4y6x2x8xy5xyΓ=f4f6f2f8f5fΓ,where y4,y6,x2,x8, and xy5, respectively, represent the field vector y in D4 and D6, x in D2 and D8, and xy in D5. A44 and A66 represent the matrix A in (2a) corresponding to the interior to interior coupling in D4 and D6. A22 and A88 represent the matrix B in (2b) corresponding to the interior to interior coupling in D2 and D8. A55 represents the matrix C in (3) corresponding to the interior to interior coupling in D5. According to the efficient WLP-FDTD method [4], A44, A66, A22, and A88 are tridiagonal matrixes. The global assembly of the interior-interface and interface-interior coupling is defined as Ai and Ai∗, respectively [5]. AΓ corresponds to the interface-interface coupling for Γ.The efficient WLP-FDTD equation (1a) is applied to the interfaces Γ1 and Γ3; therefore the coupling terms for the vector x (Γ1- the interior of D4 and Γ3- the interior of D6) do not exist in the sparse matrix system. However, the equations for vector x in D4 and D6 have the coupling term of the interior-interface Γ. According to (1b), the equation for vector x in D147 and D369 can be written as (5a)A147x147=f147-B1xyΓ1(5b)A369x369=f369-B2xyΓ3,where x147 and x369 are the vector x in D147 and D369.   D147=D1∪D4∪D7; D369=D3∪D6∪D9. xyΓ1 and xyΓ3 are the vector y in Γ1 and Γ3. B1 and B3 correspond to the interior of D147-Γ1 and the interior of D369-Γ3 coupling.According to (1a), the equations for vector y in D123 and D789 can be written as (6a)A123y123=f123(6b)A789y789=f789,where D123=D1∪D2∪D3 and D789=D7∪D8∪D9.For (4), the Schur complement system is given by [6](7)CxyΓ=g,where(8a)C=AΓ-∑iAi∗Aii-1Ai-A5∗A55-1A5i=2,4,6,8,(8b)g=fΓ-∑iAi∗Aii-1fi-A5∗A55-1f5i=2,4,6,8.Once (7) has been solved, (4) can be solved by (9a)Aiiyi=fi-AixyΓi=4,6(9b)Aiixi=fi-AixyΓi=2,8(9c)A55xy5=f5-A5xyΓ.Then (5a) and (5b) can be solved. After that, all the calculations of Eyqi,j and Exqi,j are completed.


### 2.2. Discussion about Reducing the Amount of Calculation

(5a), (5b), (6a), (6b), (9a), and (9b) are tridiagonal equations, so in this paper, the new algorithm has greatly decreased the amount of calculation compared with [2]. However, the matrix inversion is required in (8a) and (8b), and it is still very time consuming. In order to avoid solving the inverse matrix directly and further cut down the unnecessary calculation, we introduce a joint calculation scheme of (8a), (8b), (9a), (9b), and (9c). In this section, we take A4∗A44-1A4, A4∗A44-1f4, and (9a), for example.The matrixA44 can be written as(10)A44=A441A442⋱A44k⋱A44m,where m=jd-jc and n=ic-ia. A44k is a tridiagonal matrix which can be decomposed as follows:(11)A44k=LkUk=lk,1ak,2lk,2⋱⋱ak,ilk,i⋱⋱ak,nlk,n1uk,11uk,2⋱1uk,i⋱uk,n-11.According to (1a), A4 and A4∗ have only m nonzero elements. Through some complex manipulations, one can obtain(12a)A4∗A44-1A4=c1∗000000000000000,where(12b)c1∗=p1e1l1,np2e2l2,n⋱pkeklk,n⋱pmemlm,n.ek and pk represent the nonzero elements in A4 and A4∗ corresponding to A44k, respectively. (10)~(12b) indicate that the calculation of A4∗A44-1A4 is quite simple and time saving due to the efficient WLP-FDTD method in D4.Based on (10) in D4, (9a) can be written as(13)A44ky4k=f4k-A4kxyΓ1=dk.Only one nonzero elementek is contained in A4k; one can have(14)A4kxyΓ1=0,0,⋯0,ekxyΓ1,k.Sety4′k to meet the following equation:(15)A44ky4′k=f4k=d′k;then one can obtain obviously(16)A44-1f4=y4′1,y4′2,…y4′k,⋯y4′mΤ.As there is only one nonzero elementpk in A4∗ corresponding to A44k, one can have(17)A4∗A44-1f4=p1y4′1,np2y4′2,n⋯pky4′k,n⋯pmy4′m,n0⋯0Τ,where y4′k,n is the last element of y4′k in (15).Equations (13) and (15) are tridiagonal equations and the coefficient matrixes are all A44k. According to (13)~(15), only the last element is different between dk and d′k(18)dk,i′=dk,ii=1,2,…,n-1,dk,n′=dk,n+ekxyΓ1k.Therefore, we use the chasing method to the joint calculation of (13) and (15). The procedures are as follows:(19a)xk,i=dk,i-xk,i-1ak,ilk,ii=1,2,…,n-1(19b)xk,n′=dk,n′-xk,n-1ak,nlk,n(19c)y4′k,n=xk,n′(19d)xk,n=xk,n′-ekxyΓ1klk,n(19e)y4k,n=xk,n(19f)y4k,i=xk,i-uk,iy4k,i+1i=n-1,n-2,…,1.It can be seen from (19a)~(19c) that the calculation of y4′k,n is independent of xyΓ1. By inserting the results into (17), one can obtain A4∗A44-1f4. The calculation of y4k,n is dependent on xyΓ1; therefore once (7) has been solved, (19d)~(19f) can be solved. The method we introduced above avoids solving A44-1 directly. And due to the joint calculation scheme, in addition to (12a), (12b), (17), (19b), and (19c), the total amount of calculation of A4∗A44-1A4, A4∗A44-1f4, and (13) is the same as solving (13) alone. Therefore, the method has cut down the unnecessary calculation effectively. In the same way, one can deal with the other parts of ∑iAi∗Aii-1Ai, ∑iAi∗Aii-1fi, (9a), and (9b).


## 2.1. Domain Decomposition Technique

The domain decomposition technique is similar to that in [2]; the computational domain is decomposed into nine subdomains, as shown in Figure 1. However, in this paper, the conventional WLP-FDTD method is applied to D5 (fine structure exists). Meanwhile, the efficient WLP-FDTD method is used to the other subdomains and the interface Γ, where Γ=Γ1∪Γ2∪Γ3∪Γ4.Figure 1
Nine subdomains in the computational domain.In the simple and lossless medium, Maxwell’s equations of 2DTEz wave for the Laguerre-FDTD can be deduced as [4](1a)-abΔx2Eyqi-1,j+1+2abΔx2Eyqi,j-abΔx2Eyqi+1,j=2abΔxΔy∑k=0q-1Exki,j+1-Exki-1,j+1-Exki,j+Exki-1,j-2∑k=0q-1Eyki,j+2aΔx∑k=0q-1Hzki,j-Hzki-1,j-aJyqi,j(1b)-abΔy2Exqi,j-1+1+2abΔy2Exqi,j-abΔy2Exqi,j+1=-abΔxΔyEyqi+1,j-Eyqi+1,j-1-Eyqi,j+Eyqi,j-1-2∑k=0q-1Exki,j-2aΔy∑k=0q-1Hzki,j-Hzki,j-1(1c)a=2sε0,b=2sμ0,where ε0 is the free space electric permittivity, μ0 is the free space magnetic permeability, and Jyqi,j is the electric excitation source. The time-scale factor and order of the Laguerre polynomials are defined as s and q, respectively.Rewriting (1a) and (1b) as the matrix equations, we have(2a)Ayq=f1q-1(2b)Bxq=f2q-1,where vector x and y represent the unknown field variables Exqi,j and Eyqi,j, respectively.The conventional WLP-FDTD equations [1] can be written as(3)Cxyq=f3q-1,where vector xy represents the unknown field variables Eyqi,j and Exqi,j.According to (1a), in the subdomains D1~D3 and D7~D9, A in (2a) correspond to the interior to interior coupling. So firstly, yq in these subdomains can be solved by (2a) independently. Therefore one can have the sparse matrix system as(4)A440000A40A66000A600A2200A2000A880A80000A55A5A4∗A6∗A2∗A8∗A5∗AΓy4y6x2x8xy5xyΓ=f4f6f2f8f5fΓ,where y4,y6,x2,x8, and xy5, respectively, represent the field vector y in D4 and D6, x in D2 and D8, and xy in D5. A44 and A66 represent the matrix A in (2a) corresponding to the interior to interior coupling in D4 and D6. A22 and A88 represent the matrix B in (2b) corresponding to the interior to interior coupling in D2 and D8. A55 represents the matrix C in (3) corresponding to the interior to interior coupling in D5. According to the efficient WLP-FDTD method [4], A44, A66, A22, and A88 are tridiagonal matrixes. The global assembly of the interior-interface and interface-interior coupling is defined as Ai and Ai∗, respectively [5]. AΓ corresponds to the interface-interface coupling for Γ.The efficient WLP-FDTD equation (1a) is applied to the interfaces Γ1 and Γ3; therefore the coupling terms for the vector x (Γ1- the interior of D4 and Γ3- the interior of D6) do not exist in the sparse matrix system. However, the equations for vector x in D4 and D6 have the coupling term of the interior-interface Γ. According to (1b), the equation for vector x in D147 and D369 can be written as (5a)A147x147=f147-B1xyΓ1(5b)A369x369=f369-B2xyΓ3,where x147 and x369 are the vector x in D147 and D369.   D147=D1∪D4∪D7; D369=D3∪D6∪D9. xyΓ1 and xyΓ3 are the vector y in Γ1 and Γ3. B1 and B3 correspond to the interior of D147-Γ1 and the interior of D369-Γ3 coupling.According to (1a), the equations for vector y in D123 and D789 can be written as (6a)A123y123=f123(6b)A789y789=f789,where D123=D1∪D2∪D3 and D789=D7∪D8∪D9.For (4), the Schur complement system is given by [6](7)CxyΓ=g,where(8a)C=AΓ-∑iAi∗Aii-1Ai-A5∗A55-1A5i=2,4,6,8,(8b)g=fΓ-∑iAi∗Aii-1fi-A5∗A55-1f5i=2,4,6,8.Once (7) has been solved, (4) can be solved by (9a)Aiiyi=fi-AixyΓi=4,6(9b)Aiixi=fi-AixyΓi=2,8(9c)A55xy5=f5-A5xyΓ.Then (5a) and (5b) can be solved. After that, all the calculations of Eyqi,j and Exqi,j are completed.


## 2.2. Discussion about Reducing the Amount of Calculation

(5a), (5b), (6a), (6b), (9a), and (9b) are tridiagonal equations, so in this paper, the new algorithm has greatly decreased the amount of calculation compared with [2]. However, the matrix inversion is required in (8a) and (8b), and it is still very time consuming. In order to avoid solving the inverse matrix directly and further cut down the unnecessary calculation, we introduce a joint calculation scheme of (8a), (8b), (9a), (9b), and (9c). In this section, we take A4∗A44-1A4, A4∗A44-1f4, and (9a), for example.The matrixA44 can be written as(10)A44=A441A442⋱A44k⋱A44m,where m=jd-jc and n=ic-ia. A44k is a tridiagonal matrix which can be decomposed as follows:(11)A44k=LkUk=lk,1ak,2lk,2⋱⋱ak,ilk,i⋱⋱ak,nlk,n1uk,11uk,2⋱1uk,i⋱uk,n-11.According to (1a), A4 and A4∗ have only m nonzero elements. Through some complex manipulations, one can obtain(12a)A4∗A44-1A4=c1∗000000000000000,where(12b)c1∗=p1e1l1,np2e2l2,n⋱pkeklk,n⋱pmemlm,n.ek and pk represent the nonzero elements in A4 and A4∗ corresponding to A44k, respectively. (10)~(12b) indicate that the calculation of A4∗A44-1A4 is quite simple and time saving due to the efficient WLP-FDTD method in D4.Based on (10) in D4, (9a) can be written as(13)A44ky4k=f4k-A4kxyΓ1=dk.Only one nonzero elementek is contained in A4k; one can have(14)A4kxyΓ1=0,0,⋯0,ekxyΓ1,k.Sety4′k to meet the following equation:(15)A44ky4′k=f4k=d′k;then one can obtain obviously(16)A44-1f4=y4′1,y4′2,…y4′k,⋯y4′mΤ.As there is only one nonzero elementpk in A4∗ corresponding to A44k, one can have(17)A4∗A44-1f4=p1y4′1,np2y4′2,n⋯pky4′k,n⋯pmy4′m,n0⋯0Τ,where y4′k,n is the last element of y4′k in (15).Equations (13) and (15) are tridiagonal equations and the coefficient matrixes are all A44k. According to (13)~(15), only the last element is different between dk and d′k(18)dk,i′=dk,ii=1,2,…,n-1,dk,n′=dk,n+ekxyΓ1k.Therefore, we use the chasing method to the joint calculation of (13) and (15). The procedures are as follows:(19a)xk,i=dk,i-xk,i-1ak,ilk,ii=1,2,…,n-1(19b)xk,n′=dk,n′-xk,n-1ak,nlk,n(19c)y4′k,n=xk,n′(19d)xk,n=xk,n′-ekxyΓ1klk,n(19e)y4k,n=xk,n(19f)y4k,i=xk,i-uk,iy4k,i+1i=n-1,n-2,…,1.It can be seen from (19a)~(19c) that the calculation of y4′k,n is independent of xyΓ1. By inserting the results into (17), one can obtain A4∗A44-1f4. The calculation of y4k,n is dependent on xyΓ1; therefore once (7) has been solved, (19d)~(19f) can be solved. The method we introduced above avoids solving A44-1 directly. And due to the joint calculation scheme, in addition to (12a), (12b), (17), (19b), and (19c), the total amount of calculation of A4∗A44-1A4, A4∗A44-1f4, and (13) is the same as solving (13) alone. Therefore, the method has cut down the unnecessary calculation effectively. In the same way, one can deal with the other parts of ∑iAi∗Aii-1Ai, ∑iAi∗Aii-1fi, (9a), and (9b).


## 3. Numerical Results

In order to validate the effectiveness of the new algorithm, a numerical example is implemented. A parallel plate waveguide with perfect electric conductor (PEC) slot of 0.2 mm is shown in Figure2, the same as [4]. The PEC is infinitely thin. In order to simulate this fine structure, we take the graded mesh with an expansion factor α=1.07 for the y-axis (α=1.00 for the x-axis). In the x and y directions, there are 250 and 130 subdivisions including 20×Δy uniform cells for the PEC slot, and there are 10 and 40 subdivisions in D5. The size of the minimum grid is 5 mm × 10 μm. Mur’s first-order absorbing boundary conditions [7] are set to end the x boundaries. The excitation is located at the position of 80Δx,j and the wave is sinusoidal modulated Gaussian pulse:(20)Jyt=exp⁡-t-TcTd2sin⁡2πfct,where Td=1/2fc, Tc=3Td, and fc=4 GHz. The time duration Tf=4 ns, q=160, and s=3×1011 [8, 9] are chosen for the conventional WLP-FDTD and the proposed method, and Δt=1.0 ps is set to solve the Laguerre coefficient of the excitation pulse.Figure 2
Parallel plate waveguide with the PEC slot of 0.2 mm [4].Figure3 shows the vector y of the electric field at the observation point. Good agreement can be found between the conventional WLP-FDTD, the efficient WLP-FDTD [4], and the proposed method.Figure 3
Transient electric fields of the vectory at p (100Δx, 65Δy).Table1 illustrates the CPU time for the numerical simulations. It can be observed that the CPU time for the proposed method is reduced to about 15.9% of the conventional WLP-FDTD and 52.7% of the efficient WLP-FDTD. As the size of D5 is very small, the conventional WLP-FDTD calculated in this subdomain is barely time consuming. Furthermore, the iteration is not necessary in the proposed method as the splitting error does not exist in the fine structure area. Thus the efficiency has improved.Table 1
Simulation results for the different methods.
Scheme Δ t Number of iterations (q) CPU time (s) Memory (MB) WLP-FDTD 1.0 ps 161 15.57 74.58 Efficient WLP-FDTD [4] 1.0 ps 81 4.7 1.88 Proposed 1.0 ps 161 2.48 3.06Figure4 depicts the error of Ey at the observation point which is defined by Ey,error=Ey,num-Ey,WLP-FDTD, where Ey,num is the results solved by the efficient WLP-FDTD and the proposed method, respectively. Ey,WLP-FDTD is the solution utilizing the conventional WLP-FDTD. The accuracy of the proposed method is obviously improved compared with the original one.Figure 4
Error ofEy at p (100Δx, 65Δy). Calculated by the proposed method and the efficient WLP-FDTD [4].However, the memory storage required in the proposed method is slightly larger than that of the efficient WLP-FDTD. The conventional WLP-FDTD inD5 whose coefficient matrix is not a tridiagonal matrix and the Schur complement system are needful in the proposed method. So the memory storage will be increased inevitably. In this letter, all calculations have been performed on a Core 4 2.93-GHz machine.


## 4. Conclusion

In this letter, we devise a new efficient algorithm for the 2D WLP-FDTD method based on domain decomposition scheme. The whole computational domain is decomposed into several subdomains. The conventional WLP-FDTD and the efficient WLP-FDTD methods are introduced to different subdomains. And a joint calculation scheme is applied to the new algorithm. Compared with the efficient WLP-FDTD, the proposed algorithm does not need any iteration. Numerical example has verified that the proposed method can save CPU time and improve the computational accuracy. The research for the 3D case will be presented in the further work.



---

*Source: 3163781-2016-11-30.xml*