Domain Decomposition CN-FDTD with Unsplit-Field PML for Time-Reversed Channel Analysis

In this paper, an efficient domain decomposition (DD) technique is introduced into the implicit CrankNicolson finite-difference time-domain (CN-FDTD) method to analyze the channel characteristics of time reversal (TR) waves. As an unconditionally stable time-marching method, DD-CN-FDTD is suitable to solve the multiscale problem involving special microstructures in a sub-wavelength array. The standard unsplit-field perfectly matched layer, which is implemented with auxiliary differential equations, is derived here to truncate the computational domain of a multipath indoor environment. Furthermore, each sub-matrix is preconditioned by the reverse CuthillMckee scheme for easier lower-upper decomposition. Numerical results of TR wave propagation demonstrate the performance of the proposed method.


Introduction
The finite-difference time-domain (FDTD) method has been widely used for solving electromagnetic problems, but its time step is constrained by the Courant-Friedrich-Levy (CFL) stability condition [1].To remove this limit, the Crank-Nicolson (CN) scheme is introduced into FDTD to achieve unconditional stability [2], and it is more efficient than the conventional FDTD for solving multiscale problems.Recently, a domain decomposition (DD) technique has been introduced into weighted Laguerre polynomials (WLPs) based FDTD to get efficient solutions [3], [4].
The ultrawide-band (UWB) communication has gathered great interest due to its improved multiple-access capabilities and immunity to interference, but the communication quality will be decreased due to reflection, diffraction and scattering in a multipath propagation environment.
Fortunately, the time reversal (TR) signal processing technique in [5] can make full use of multipath propagation created by various scatterers in an indoor environment [6].In the TR scheme, the received power is concentrated at a specific time and location and very low co-channel interference in a multipath system can be achieved due to spatial and temporal focusing [7], [8].For the numerical simulation of TR wave propagation, the number of unknowns is often large due to the fine structures and big computational range.Therefore, the generated large-scale sparse matrix in CN-FDTD leads to a heavy computational burden.
In this paper, the DD-CN-FDTD method is employed to simulate the multiscale electromagnetic problem of TR channel analysis in a two-dimensional (2-D) multipath indoor environment.The standard unsplit-field perfectly matched layer (PML) is derived for the absorbing boundary condition of DD-CN-FDTD with auxiliary differential equations (ADEs) [9].Moreover, the reverse Cuthill-Mckee (RCM) technique [10] is employed to compress the bandwidth of sparse matrices and largely reduce the complexity of lower-upper (LU) decomposition.A numerical example of TR wave propagation is calculated to show the potential application of the proposed method.
Using the central difference scheme, we can transform the differential equations in (1) into difference ones.Then inserting (1c) into (1a) and (1b) with reference to [2], respectively, we can obtain the implicit updating equations for the electric field in CN-FDTD.Rewriting the implicit equations of E x and E y as a matrix form, we get where A is the coefficient matrix, E x,y is the unknown field vectors, b is the known terms and J is the current density.

Implementation of Domain Decomposition
In the DD structure, for example, the whole computational domain is decomposed into two subdomains, with the interface between the subdomains named as Γ 1_2 , as shown in Fig. 1.For a specific simulated structure, a nonconformal DD scheme in which different sub-domains are discretized by different grids has been employed in the weighted Laguerre polynomial based FDTD [11].The nonconformal DD scheme can also be implemented in CN-FDTD to handle the different grids on the interface between neighboring subdomains.Here, for simplicity, the graded nonuniform grids are used to divide the whole computational domain and the conformal domain interfaces are used.
The unknowns, namely the field components on the grids, are reordered by starting with the ones in D 1 , followed by those in D 2 , and ending with those associated with Γ 1_2 .Then, the large matrix system (2) generated by the 2-D CN-FDTD method can be rewritten as [4] Without loss of generality, it is assumed that the whole computational domain is decomposed into N subdomains.The matrix equation for the N-subdomain solution has the same form as (3).From (3), the linear system that is local to the ith subdomain and the linear system that is related to the interface can be written as , , 1 From both (4) and ( 5), the variable E x,yi can be eliminated from ( 5) and the Schur complement can be got where Once the interface variables are obtained by solving the Schur complement system ( 6), ( 3) can be decoupled by solving the following subdomain problems , .
It is evident that the systems of ( 7) are independent of each other and can be solved in a parallel manner.

Formulations of Unsplit-Field PML
The frequency domain Maxwell's equations of a 2-D TE z wave in complex stretched coordinates can be written as [9] , , 0 , , , , 0 , , where j 1   .And S v (v = x, y) are the stretched coordinate matrix coefficients defined via the complex frequency shifted (CFS) PML parameters 0 , , +j where α v , κ v and σ v are assumed to be positive real, and they are one-dimensional functions along the v-axis.
Particularly, when α v = 0 and κ v = 1, the standard unsplitfield PML can be derived for CN-FDTD.
To simplify the problem, but without loss of generality, the E x component is taken as an example.With reference to [9], (8a) can be expressed as , , , , The auxiliary variable is introduced into (10) as , , , , , 0 and then the expression of ( 10) can be directly transformed into a time-domain form , , , , , 0 , , where the auxiliary variable Similar to the derivation of ( 12) and ( 13), the timedomain expressions of other electromagnetic field components and auxiliary variables can also be obtained.
With reference to [2], the 2-D implicit E x formulation of the standard unsplit-field PML for CN-FDTD can be given by using the central difference scheme to get the difference forms of ( 12) and (13).Inserting (13) and the difference equation of H z into (12), the implicit updating equation for E x can be got.Similarly, the updating equation of E y can also be obtained.

Bandwidth Compression by RCM
Since the banded-sparse coefficient sub-matrices produced by DD-CN-FDTD keep unchanged during the timemarching process, the LU decomposition for each submatrix needs to be performed only once at the beginning of the calculation.The computational complexity of LU decomposition is O(2m 1 m 2 n + m 2 n -2n), where n is the number of unknowns, m 1 is the bandwidth of the lower triangular matrix, and m 2 is the bandwidth of the upper one.The memory requirement of LU decomposition is O(kn), where k is the bandwidth of the coefficient matrix [12].In order to save calculation time and reduce memory requirement, the RCM technique [10] that generates a symmetric sparsity pattern with a small bandwidth is used for the coefficient matrices generated from the sub-domains.

Numerical Results
To validate the accuracy and efficiency of the proposed DD-CN-FDTD method for solving multiscale and electrically large electromagnetic problems, a numerical example of TR wave propagation is calculated in this section.In this section, a 2-D structure is analyzed for simplicity.However, the proposed method in this paper can be extended to solve three-dimensional scenarios.Here, all calculations are performed on an Intel (R) Core (TM) i3-2120×4 3.30 GHz machine with 12 GB RAM.

Channel Characteristics Extraction
To begin with, the channel characteristic of TR waves in a complicated multipath indoor environment is considered, as shown in Fig. 2.There are several fixed perfectly electric conductor (PEC) objects and one moving PEC object in the computational model.The box and chairs are dielectrics with ε r = 1.5, and the walls are with ε r = 4.A sub-wavelength array consists of the special microstructure of several metal sheets and it is loaded to the excitation source.The thickness of each sheet is 0.5 mm and its length is 400 mm.
To accurately model the geometry of the microstructure, the graded nonuniform grids along the yaxis are employed to divide the whole computational domain.The 2-D domain is 6 m × 4 m and is meshed into 300 × 280 rectangular grids, with the minimum cell size of 20 mm × 0.25 mm and the maximum one of 20 mm × 20 mm.The whole domain is terminated by 10 cells thick unsplit-field PML.The time variation of the source is given by where f 0 = 0.65 GHz, f max = 1 GHz, τ = 1 / (2f max ), and t 0 = 3τ.The time step for FDTD is assigned to be Δt FDTD = 0.4167 ps which satisfies the CFL condition, while the time steps for CN-FDTD and DD-CN-FDTD are Δt CN-FDTD = Δt DD-CN-FDTD = 100 Δt FDTD .In the DD structure, the computational domain is decomposed into four subdomains.
For the first procedure of TR simulation, the excitation is injected into H z of the Source node.Conventional FDTD consists of 100 000 time steps, while both DD-CN-FDTD and CN-FDTD consist of 1 000 time steps, covering the same time interval.For the multiscale problems, it is computationally challenging for the conventional FDTD method due to the CFL stability limit.However, DD-CN-FDTD is an unconditionally stable scheme, and its time step is no longer restricted by the CFL stability limit.Thus, the much larger time step of DD-CN-FDTD results in the much smaller time-marching number.The UWB signal of H z is recorded at Probe node.Then, the channel impulse response (CIR) of received UWB signal is extracted from the time-domain waveform by using the CLEAN algorithm presented in [13], as shown in Fig. 3.In a multipath indoor environment, the discrete impulse response h(t) of the channel is modeled as a summation of reflected components.The real reflection response y(t) to the input signal s in (t) can be computed by convolving s in (t) with h(t).Then, h(t) of a UWB channel can be extracted by deconvolving the reflection signal y(t) from s in (t) with the CLEAN algorithm.
In the second procedure of TR simulation, the timereversed waveform of the UWB signal, serving as the excitation function, is injected from the Probe node back into the whole space.N t = 112 000 for FDTD while N t = 1 120 for DD-CN-FDTD and CN-FDTD are assigned.
Here, the time-marching number should be larger than that in the first procedure to ensure the time-reversed waveform of the UWB signal recorded in the first procedure is totally injected into the whole space.The space-time focusing can be realized at the original Source node and the TR-UWB signal is recorded.The CIR of TR-UWB signal is also extracted by CLEAN, as shown in Fig. 4. From Figs. 3 and  4, the results from two methods are in good agreement.
To quantify the UWB and TR-UWB channels, several parameters of channel characteristics are defined to evaluate the performance of TR in wireless communication [14], [15].A smaller root mean square (RMS) delay spread δ means reduced intersymbol interference (ISI) and improved communication quality, and a larger signal to noise ratio (SNR) R leads to better signal quality.Furthermore, the positive value of the temporal compression parameter C RMS indicates that the RMS delay spread of the channel is reduced by TR.If the gain of signal to noise ratio G SNR is larger than 1, the SNR of the channel is also improved by the TR technique.Table 1 shows the UWB channel characteristic parameters calculated from the three methods, and Table 2 shows the TR-UWB ones.It can be revealed from Tabs. 1 and 2 that the results from DD-CN-FDTD agree quite well with those from FDTD and CN-FDTD.As demonstrated in Tab. 2, C RMS is positive and G SNR is larger than 1, indicating that the TR technique reduces the ISI and improves the SNR due to its time-space focusing effect.
A comparison of computing resources among the three methods is shown in Tab. 3. DD-CN-FDTD results in a faster calculation with the pre-conditional RCM technique than the other two methods.Furthermore, the memory requirement of matrices L and U in DD-CN-FDTD is smaller than that in CN-FDTD.

Time-varying Communication System
The performance of the TR technique in a time-varying channel environment with fixed terminals (Source and Probe) and fixed objects is performed.A PEC bulk is mov-

Conclusion
In this paper, the standard unsplit-field PML with ADEs is proposed for CN-FDTD and a DD technique is introduced to efficiently solve multiscale and electrically large TR problems.Compared with CN-FDTD, the smallscale sub-matrices in DD-CN-FDTD result in fast calculation and low memory requirement.The proposed method is used to calculate channel characteristic parameters of TR waves in a multipath indoor environment, and the results demonstrate its accuracy and efficiency.