International Journal of Antennas and Propagation
Volume 2019 (2019), Article ID 5043028, 11 pages
https://doi.org/10.1155/2019/5043028

Research Article

Analysis of GPR Wave Propagation in Complex Underground Structures Using CUDA-Implemented Conformal FDTD Method

Jianwei Lei,¹ Zibin Wang,2,3 Hongyuan Fang,¹ Xin Ding,3,4 Xiaowang Zhang,⁵ Man Yang,¹ and Haitao Wang3,6

¹College of Water Conservancy & Environmental Engineering, Zhengzhou University, Zhengzhou 450001, China
²Guangdong Hualu Communications Technology Co., Ltd., Guangzhou 510420, China
³National Local Joint Engineering Laboratory of Major Infrastructure Testing and Rehabilitation Technology, Zhengzhou 450001, China
⁴Airport Management Office of Henan Toll Expressway Management Center, Zhengzhou 450019, China
⁵Henan Transportation Research Institute Co., Ltd., Zhengzhou 450006, China
⁶Zhengzhou Development & Investment Group Co., Ltd., Zhengzhou 450000, China

Correspondence should be addressed to Hongyuan Fang; 18337192244@163.com

Received 29 March 2019; Revised 4 June 2019; Accepted 17 June 2019; Published 19 August 2019

Academic Editor: Luciano Tarricone

Copyright © 2019 Jianwei Lei et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Ground penetrating radar (GPR), as a kind of fast, effective, and nondestructive tool, has been widely applied to nondestructive testing of road quality. The finite-difference time-domain method (FDTD) is the common numerical method studying the GPR wave propagation law in layered structure. However, the numerical accuracy and computational efficiency are not high because of the Courant-Friedrichs-Lewy (CFL) stability condition. In order to improve the accuracy and efficiency of FDTD simulation model, a parallel conformal FDTD algorithm based on graphics processor unit (GPU) acceleration technology and surface conformal technique was developed. The numerical simulation results showed that CUDA-implemented conformal FDTD method could greatly reduce computational time and the pseudo-waves generated by the ladder approximation. And the efficiency and accuracy of the proposed method are higher than the traditional FDTD method in simulating GPR wave propagation in two-dimensional (2D) complex underground structures.

1. Introduction

GPR is a nondestructive way to detect the internal material distribution of underground structure using electromagnetic waves [1]. Its advantages include being fast, continuous, high resolution and no damage. Because of these advantages, it has been widely applied to geological prospecting, civil engineering, highway reconnaissance, and many other areas [2–6]. Numerical simulation of GPR electromagnetic wave propagation in underground structures is helpful to better explain the real GPR profiles. It can provide a way to explore the link between subsurface structures properties and GPR data.

Various approaches have been developed to simulate GPR wave propagation. These include the ray-tracing method [7], the transmission-line-matrix method [8], the finite element method (FEM) [9], the semianalytic mode matching algorithm [10], the pseudo-spectral method [11], the FDTD method [12], the Alternating Direction Implicit- (ADI-) FDTD method [13], and the symplectic method [14]. Among these, the FDTD method proposed by Yee in 1966 [15] has been the most popular simulation tool due to its great stability and convergence. It uses time-step finite-difference to approximate the differential form of Maxwell’s equations for Faraday and Ampere’s laws and calculates the magnetic and electric field vector components at regularly arranged, distinct points in space. The current state of the art in GPR-oriented FDTD approaches has developed rapidly; for example, a software tool called GprMax has been successfully employed in the detection mechanism of GPR [16], a GPU-accelerated FDTD solver integrating it into open source EM simulation software was developed to modeling GPR [17], and a systematic framework for accurate and realistic numerical modeling of GPR for landmine detection has been developed [18]. However, the FDTD method requires large memory and calculation time in simulating GPR wave propagation in 2D complex underground structures due to the limitations of the CFL stability condition [19].

The accuracy and efficiency of forward model are not only related to the performance of the algorithm itself, but also related to computer acceleration technology and modeling method. There are three main methods for improving the accuracy when determining the electromagnetic response of a simulated irregular shape target: one is to increase the density of the discrete grids; however the computation and time also increase correspondingly. For complex problems, the computational resources are difficult to meet. Another approach is to adopt subgrid technology, such that the irregular shape target area uses a subdivision grid, and a normal grid is adopted for other areas. The subgrid technology can greatly improve the calculation accuracy and efficiency, but it is also easy to generate reflection waves at the interface of the coarse and fine grids. A last approach is conformal grid technology [20]. This approach can deal with objects with arbitrary surfaces, and the formula deductions and calculations are simple.

In 2006, NVIDIA developed a general-purpose parallel computing architecture called compute unified device architecture (CUDA) that leverages parallel compute engine in NVIDIA graphic processor unit (GPU) to solve many complex computational problems [21]. Due to its excellent programmability and remarkable computational power for parallel executions, CUDA has gained tremendous popularity in the scientific computing society. GPU computing has developed in various fields, such as acoustics [22], biomedical applications [23], plasmonics [24], dispersive media [25], and computational electromagnetism [26–34]. FDTD based on the GPU acceleration technique has been applied in the GPR simulation model [35, 36], but the pseudo-waves that are generated by the ladder approximation in FDTD modeling method are not considered. This study combines FDTD method and surface conformal technique and GPU acceleration technique proposing a precise and efficient forward modeling algorithm of GPR which can greatly reduce computation time and the pseudo-waves that are generated by the ladder approximation.

In this paper, the simulation model of GPR wave propagation in 2D underground structure was established by CUDA-implemented conformal FDTD method. Line source was chosen to be the incident wave, and second-order Mur absorbing boundary condition was adopted to avoid reflections from the modeled edges. Three numerical examples were provided to verify the high efficiency and accuracy of the CUDA-implemented conformal FDTD method. One obtained the simulated reflection waveform from the two-layer pavement model, the other got the simulated GPR trace profile of the pavement model with structural damage, and the last one demonstrated the GPR simulation section of underground structures with voids.

2. FDTD Method

2.1. The Governing Equations

Among Maxwell’s electromagnetic equations, the FDTD method uses two curl equations:

Here is the electric field, is the magnetic field, ε is the dielectric constant, σ is the conductivity, and μ is the permeability.

Considering 2D Transverse Magnetic (TM) case, (1) can be written as

For 2D TM case, the distribution of the electric and magnetic fields of Yee’s cell of the 2D TM wave is shown in Figure 1 [37].

Figure 1: Schematic diagram of position of the electric and magnetic field vector of the Yee TM wave.

By applying the central finite-difference approximation to the spatial and time derivatives of (2) yields, the 2D TM wave discrete equations can be written as

where , .

2.2. Mur Absorbing Boundary Conditions

An absorbing boundary condition allows the FDTD to calculate approximately the infinite domain problems with limited computer capacity. Its basic idea is to truncate boundaries in the computing area through special treatment or add the ideal mathematical absorption material to create the incident wave on the interface of the electromagnetic wave reflection phenomenon. The absorption boundary applied in this paper is the second-order Mur absorption boundary [38].

For the TM wave, the distribution map of the field component at the left truncated boundary is shown in Figure 2, and the second-order Mur absorption condition on the left edge can be expressed as

Figure 2: The arrangement of field components at the left truncating boundary.

Equation (3) is discrete at point . Using a partial differential in the central difference substitution formula, the resulting difference scheme can be expressed as

The points on other boundaries and corners can be calculated using the same method.

2.3. The Stability Condition

Numerical stability is enforced by the Courant condition, which limits the ratio between the spatial discrete step and time step in the simulation. In 2D grid, this condition is given in the following inequality:

where c is the speed of light, Δt is the time step, and Δx and Δy are the spatial step in the x and y direction, respectively.

2.4. Conformal Grid Technology

We adopted a conformal grid method based on effective medium parameters to simulate an arbitrary curved boundary of subsurface structure. A two-dimensional circle is shown as an example to explain the conformal grid method theory in Figure 3. Figure 3(a) shows the actual subdivision grid diagram for a circle, the white grids are normal grid points and the green grid is actual subdivision grid of a circle. Figure 3(b) shows the subdivision grid diagram for the conventional ladder approximation of a circle, and the blue grid is conventional ladder approximation grid of a circle. Figure 3(c) shows the conformal grid subdivision diagram for the circle, and the purple grid is conformal grid. The grid points that are entirely inside the circle or entirely outside of the circle are classified as normal grid points, while the rest are classified as conformal grid points. The effective medium parameters are adopted, to process the conformal grids in this paper.

(a)

(b)

(c)

Figure 3: Schematic view for the subdivisions of a circle. (a) An actual model for a circle. (b) The subdivision diagram using the conventional ladder approximation. (c) The conformal grid for a circle.

Next, a conformal cell is extracted from Figure 3(c) to illustrate the effective medium parameters for the conformal grid points in the two-dimensional TM wave example. In the FDTD methods, the electric field E is located at the midpoint of the conformal cell, and the magnetic field H is located at the edge of the conformal cell. The calculation for a conformal grid point of the medium in Cartesian coordinates is shown in Figure 4, where N is the sampling point of the field components ; A and B are the sampling points of the field components ; C and D are the electric field sample points of the field components . Δx and Δy are the width and height of a grid, and and are the length of the different media on a corresponding edge; and are the area of media 1 and 2, respectively. We assume that the electromagnetic parameters of media 1 and 2 are ε₁, σ₁, μ₁ and ε₂, σ₂, μ₂, respectively.

Figure 4: Schematic view for the equivalence of the medium parameters for a conformal grid point. (a) An actual subdivision grid. (b) A conformal grid.

As the magnetic field node is located at the middle point of the grid edge, the effective values of the magnetic conductivity and permeability at each magnetic field node can be calculated by a weighted average of the length of the different media on a corresponding edge. Similarly, the effective values of the dielectric constant and conductivity at each electric field node can be calculated by a weighted average of the area occupied by different media on the corresponding surface. The effective magnetic permeability of the A, B, C, and D nodes in Figure 4(b) can be expressed as

The effective magnetic permeability of the conformal mesh, respectively, can be written as

The effective dielectric constant and the effective conductivity of the N point in Figure 4(b) can be shown as

where eff denotes the equivalent value of the parameter.

Substituting (7)-(9) into (3), we obtain the differential iterative equations for the conformal FDTD method:

where , .

3. The Simulation of GPR Waves on GPU

3.1. CUDA-Implemented FDTD Method

CUDA programming is well suited to solve problems that can be expressed as parallel computing, which achieves application acceleration performance by mapping data elements to multiple threads. The kernel is a function that is executed on the device, which is the most important component of the CUDA program. The kernel is always represented as a sequential program on the host and executed on the device. CUDA exposes a two-level thread hierarchy that provides an optimization method for developers. All threads are divided into blocks, all of which are grouped into a single grid. The kernel executes in the grid of thread blocks.

E and H field components calculate in a bidimensional xy square domain. The number of threads is decided by the calculation domain, which is the number of Yee units, because we map threads to units. As shown in Figure 5, the strategy adopted is mapping each thread to Yee’s cell and each thread block to a group of cells. Each thread deals with one Yee cell during the parallel computing. In this way, the entire calculation domain is wrapped by the CUDA grid.

Figure 5: Mapping between real domain and CUDA grid of thread blocks.

In 2D case, the flowchart of CUDA kernel performing the field calculation is depicted in Figure 6 [39]. There are two CUDA kernels performing the field calculations. The first kernel computes E, and the second updates H. T is time step; is total time step.

Figure 6: The flowchart of the CUDA-implemented FDTD method.

3.2. The Forward Modeling

The steps of the forward modeling of the propagation of GPR in the layer pavement system are as follows:

(1)Input basic parameters of pavement structure and build pavement structure model(2)Obtain and input the source of excitation; make numerical difference a discrete source and make it a point corresponding to the time step(3)Calculate the propagation characteristics of GPR electromagnetic waves in the pavement structure with the finite-difference principle of FDTD, including reflection, loss characteristics, analysis, and calculation of the general field(4)Calculate the reflection field(5)Output the signal of radar reflection

According to the above steps, a forward model of the GPR propagation in the pavement structure is established. The block diagram of the program is shown in Figure 7.

Figure 7: The forward modeling procedure of GPR electromagnetic wave in the pavement structures.

3.3. The Accuracy of the Proposed Algorithm

The following numerical simulations were performed on a desktop computer with Intel Core i7-6700K CPU and NVIDIA Geforce GTX 1070 Graphics card. The desktop operating system is Windows 7 Pro 64-bit. The development tool of CUDA program is the CUDA Toolkit 7.5 built with Microsoft Visual Studio 2010. In the following, we verified the accuracy of the absorption boundary and algorithm by a simple example.

In the 2D FDTD region in free space, the source of excitation was a sine point source E(t) = sin(2π · t · dt · m) (m= 1.5e13), and Mur boundary was used at the regional truncation boundary. In order to ensure the stability and convergence of FDTD method, the spatial discrete grid and the time step were chosen as Δx = Δy = 0.25 cm and dt = 0.005 ns, respectively. The point source was initiated in the center of the 50 cm × 50 cm square free space by using the 2D TM wave as an example.

A CUDA calculation model was established for each thread to solve a Yee cellular (magnetic) field. There were 16 × 16 threads in each block, and the simulated area was divided into 200 × 200 grids. The point source and observation point were located at (25 cm, 25 cm) and (25 cm, 30 cm), respectively. Using traditional FDTD method and CUDA-implemented FDTD method, we calculated the 600-time-step field change with the same parameter calculation. As shown in Figure 8, the simulation results of traditional FDTD method and CUDA-implemented FDTD method were with excellent agreement.

Figure 8: Field value of the observation point with two methods.

As shown in Figure 9, it was the 600-time-step field change diagram calculated by the CUDA-implemented FDTD method and traditional FDTD method. It can be seen that GPR electromagnetic wave reached the boundary without causing reflection due to the setting of absorbing boundary. It also simulated intuitively the GPR electromagnetic wave propagation in free space and proved the accuracy of the Mur absorbing boundary.

(a)

(b)

Figure 9: Variation diagram of the field value. (a) Traditional FDTD method. (b) CUDA-implemented FDTD method.

4. Numerical Simulation

Firstly, we verified the high efficiency of the CUDA-implemented FDTD method by simulating a two-layer pavement model. As shown in Figure 10, the top layer of the model represented air with = 1, σ = 0; the second layer represented asphalt concrete with = 6, σ = 1 mS / m; the third layer was the base with = 30 and σ = 2 mS / m. We assumed that all materials are nonmagnetic. The time step and spatial increment were set equal to dt = 0.005 ns and Δx = Δy = dl = 0.25 cm, respectively. We used a Ricker wavelet with a unit amplitude and a center frequency of 1.0 GHz (Figure 11) to represent the line source of the GPR transmitter.

Figure 10: Computational domain of the three-layered model.

Figure 11: Map of Ricker waveform.

The transmitter was at the point (15 cm, 100 cm) and the receiver was at the point (15 cm, 110 cm), respectively. The simulation results of traditional FDTD method and CUDA-implemented FDTD method after 1000 time steps were shown in Figure 12. Excellent agreement was reached. It was proved that traditional FDTD method and CUDA-implemented FDTD method have the same accuracy. Moreover, the CUDA-implemented FDTD method consumed 26.276 s, while the traditional FDTD method increased to 299.919 s. The former saved 91.24% of the time. The result demonstrated that CUDA-implemented FDTD method greatly reduced computation time.

Figure 12: Map of the single-track radar data calculated using two methods.

As shown in Figure 12, it can prove that CUDA-implemented FDTD method has the same accuracy as the traditional FDTD method. The speedup ratios achieved are shown in Table 1. When the number of units is small, the acceleration ratio is also small, so the acceleration performance obtained is not significant. When the number of units is getting larger, better acceleration performance will be obtained. In particular, when the FDTD cells are up to 4096×4096, we can obtain more than 13 times speedup ratio.

Table 1: Speedup ratio of the two methods.

FDTD Cells
CPU Times (s)
GPU Times (s)
Speedup Ratio

128 × 128
1.063
0.398
2.67

256 × 256
4.290
0.743
5.77

512 × 512
18.981
1.879
10.11

1024 × 1024
72.834
6.773
10.75

1536 × 1536
158.545
15.016
10.56

2048 × 2048
299.919
26.276
11.41

2560 ×2560
475.560
41.643
11.58

3072 × 3072
701.439
59.439
11.81

4096 × 4096
1357.060
105.684
13.02

GPR profiles of a structurally damaged pavement were simulated in second example. As shown in Figure 13, the first layer of the model was an air layer, and the bottom was three-layer pavement with structural damage. There was a 10 mm wide crack in first layer. There were a circular no-dense region with radius of 5 cm and a 30 cm × 40 cm rectangular cavity in second layer. In this example, assuming that the porosity of the region was 10%, all dielectric parameters of the material were shown in Figure 13. Receivers and transmitters were placed along the air-earth interface every 1 cm to measure the reflection waveform of the GPR wave in the pavement structure. The distance between receiver and transmitter was 5 cm. The incident wave in first example was used to represent the line source of the transmitter in this model. The increments of the time step and spatial interval were selected as dt = 0.01 ns and Δx = Δy = dl = 0.5 cm, respectively. GPR tracking profiles of the model simulated by CUDA-implemented FDTD method and traditional FDTD method were plotted in Figures 14 and 15.

Figure 13: Schematic diagram of a two-layer pavement model with no-dense region.

Figure 14: GPR trace profile with traditional FDTD method.

Figure 15: GPR trace profile with CUDA-implemented FDTD method.

As shown in Figures 14 and 15, a good agreement has been reached. However, the traditional FDTD method required 2.7357 h; the CUDA-implemented FDTD method required 0.1771 h. The results showed that the reflection waveform of CUDA-implemented FDTD method was consistent with the traditional FDTD method, but the CUDA-implemented method reduced the computational time by more than 93.5%.

In this example, a complex underground structure with two cavities was established. We verified the high efficiency and accuracy of CUDA-implemented conformal FDTD method by using conformal and nonconformal modeling method for two cavities, respectively. As shown in Figure 16, this model included an air-earth interface and a soil layer with two cavities with a radius of 5 cm. The model size and dielectric parameters were shown in Figure 16. The excitation source was the same as the second model. The time step and spatial interval increment were selected as dt = 0.01 ns and Δx = Δy = dl = 0.5 cm, respectively.

Figure 16: Sketch of the complex underground structure model.

We placed receivers and transmitters along the air-earth interface every 1 cm. The distance between the receiver and transmitter was 10 cm. Figures 17 and 18 were GPR simulation section using CUDA-implemented FDTD method and traditional FDTD method. The right cavity used conformal modeling method, the left one with nonconformal modeling method. And from Figures 17 and 18, it can be seen that the conformal modeling method has fewer pseudo-waves than the nonconformal modeling method. The CUDA-implemented FDTD method required 967.23 s, and traditional FDTD method required 11606.76 s. The results indicate that CUDA-implemented conformal FDTD algorithm can greatly reduce computational time and the pseudo-waves that are generated by the ladder approximation.

Figure 17: GPR simulation section with CUDA-implemented FDTD method.

Figure 18: GPR simulation section with traditional FDTD method.

5. Conclusions

This paper presented a CUDA-implemented conformal FDTD method for simulating GPR wave propagation in 2D low-loss medium. The high efficiency and accuracy of the method were verified by three numerical examples. The results showed that CUDA-implemented conformal FDTD algorithm can greatly reduce computation time and the pseudo-waves that were generated by the ladder approximation. In addition, we obtained the GPR simulation profile of complex underground structure by forward modeling using the CUDA-implemented conformal FDTD method. By analyzing the simulated profiles of GPR, we can deepen our understanding of the propagation of GPR waves in underground space and better interpret the real GPR profiles. This paper focuses on 2D isotropic conditions. However, this method can be extended to three-dimensional (3D) or anisotropic situations. Through further research, we will develop a CUDA-implemented conformal FDTD to simulate the propagation of GPR waves in anisotropic media or 3D structures.

Data Availability

The data used to support the findings of this study are included within the article.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

This research was supported by the National Natural Science Foundation of China (No. 51678536, No. 41404096), the National Key Research and Development Program of China (No. 2017YFC1501200), the Scientific and Technological Research Program of Henan Province (No. 152102310066), and the Henan Transportation Science and Technology Project of Henan Province in 2018 (No. 2018J7) for which the authors are grateful.

References

Q.-W. Dai, D.-S. Feng, and J.-S. He, “Finite difference time domain method forward simulation of complex geoelectricity ground penetrating radar model,” Journal of Central South University of Technology (English Edition), vol. 12, no. 4, pp. 478–482, 2005.

N. J. Cassidy and T. M. Millington, “The application of finite-difference time-domain modelling for the assessment of GPR in magnetically lossy materials,” Journal of Applied Geophysics, vol. 67, no. 4, pp. 296–308, 2009.

U. Böniger and J. Tronicke, “Integrated data analysis at an archaeological site: A case study using 3D GPR, magnetic, and high-resolution topographic data,” Geophysics, vol. 75, no. 4, pp. B169–B176, 2010.

Z. Zyada, T. Matsuno, Y. Hasegawa, S. Sato, and T. Fukuda, “Advances in GPR-based landmine automatic detection,” Journal of The Franklin Institute, vol. 348, no. 1, pp. 66–78, 2011.

C. Anitescu, Y. Jia, Y. J. Zhang, and T. Rabczuk, “An isogeometric collocation method using superconvergent points,” Computer Methods Applied Mechanics and Engineering, vol. 284, pp. 1073–1097, 2015.

K. M. Hamdia, H. Ghasemi, X. Zhuang, N. Alajlan, and T. Rabczuk, “Sensitivity and uncertainty analysis for flexoelectric nanostructures,” Computer Methods Applied Mechanics and Engineering, vol. 337, pp. 95–109, 2018.

X. Zeng, G. A. McMechan, J. Cai, and H. W. Chen, “Comparison of ray and Fourier methods for modeling monostatic ground- penetrating radar profiles,” Geophysics, vol. 60, no. 6, pp. 1727–1734, 1995.

C. Liu and L. Shen, “Numerical simulation of subsurface radar for detecting buried pipes,” IEEE Transactions on Geoscience and Remote Sensing, vol. 29, no. 5, pp. 795–798.

A. Villarino, B. Riveiro, D. Gonzalez-Aguilera, and L. J. Sánchez-Aparicio, “The integration of geotechnologies in the evaluation of a wine cellar structure through the finite element method,” Remote Sensing, vol. 6, no. 11, pp. 11107–11126, 2014.

A. W. Morgenthaler and C. M. Rappaport, “GPR wave scattering from complex objects using the semi-analytic mode matching algorithm: coordinate scattering center selection,” IEEE Transactions on Geoscience and Remote Sensing, vol. 49, no. 6, pp. 1949–1956, 2011.

J. M. Carcione, S. Valle, and G. Lenzi, “GPR modelling by the Fourier method: improvement of the algorithm,” Geophysical Prospecting, vol. 47, no. 6, pp. 1015–1029, 1999.

E. Karpat, M. Çakır, and L. Sevgi, “Subsurface imaging, FDTD-based simulations and alternative scan/processing approaches,” Microwave and Optical Technology Letters, vol. 51, no. 4, pp. 1070–1075, 2009.

D.-S. Feng and Q.-W. Dai, “GPR numerical simulation of full wave field based on UPML boundary condition of ADI-FDTD,” NDT & E International, vol. 44, no. 6, pp. 495–504, 2011.

T. Hirono, W. W. Lui, and K. Yokoyama, “Time-domain simulation of electromagnetic field using a symplectic integrator,” IEEE Microwave and Guided Wave Letters, vol. 7, no. 9, pp. 279–281, 1997.

K. S. Yee, “Numerical solution of initial boundary value problems involving Maxwell equations in isotropic meidia,” IEEE Transactions on Antennas and Propagation, vol. 14, no. 3, pp. 302–307, 1966.

A. Giannopoulos, “Modelling ground penetrating radar by GprMax,” Construction and Building Materials, vol. 19, no. 10, pp. 755–762, 2005.

C. Warren, A. Giannopoulos, A. Gray et al., “A CUDA-based GPU engine for gprMax: Open source FDTD electromagnetic simulation software,” Computer Physics Communications, vol. 237, pp. 208–218, 2019.

I. Giannakis, A. Giannopoulos, and C. Warren, “A realistic FDTD numerical modeling framework of ground penetrating radar for landmine detection,” IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, vol. 9, no. 1, pp. 37–51, 2016.

A. Taflove, Computational Electrodynamics: The Finite-Difference Time-Domain Method, Artech House, Norwood, MA, USA, 3rd edition, 2005.

X.-J. Hu, D.-B. Ge, and B. Wei, “Study on MCFDTD for 3-D coated targets by using effective parameters,” Journal of Systems Engineering and Electronics, vol. 28, no. 11, pp. 1652–1667, 2006.

E. Lindholm, J. Nickolls, S. Oberman, and J. Montrym, “NVIDIA Tesla: a unified graphics and computing architecture,” IEEE Micro, vol. 28, no. 2, pp. 39–55, 2008.

J. J. Lopez, D. Carnicero, N. Ferrando, and J. Escolano, “Parallelization of the finite-difference time-domain method for room acoustics modelling based on CUDA,” Mathematical and Computer Modelling, vol. 57, no. 7, pp. 1822–1831, 2013.

. Jieru Chi, . Feng Liu, E. Weber, . Yu Li, and S. Crozier, “GPU-accelerated FDTD modeling of radio-frequency field–tissue interactions in high-field MRI,” IEEE Transactions on Biomedical Engineering, vol. 58, no. 6, pp. 1789–1796, 2011.

K. H. Lee, I. Ahmed, R. S. Goh, E. H. Khoo, E. P. Li, and T. G. Hung, “Implementation of the fdtd method based on lorentz-drude dispersive model on gpu for plasmonics applications,” Progress in Electromagnetics Research, vol. 116, pp. 441–456, 2011.

M. R. Zunoubi, J. Payne, and W. P. Roach, “CUDA implementation of TEZ-FDTD solution of Maxwell's equations in dispersive media,” IEEE Antennas and Wireless Propagation Letters, vol. 9, pp. 756–759, 2010.

N. Takada, N. Masuda, T. Tanaka, Y. Abe, and T. Ito, “A GPU implementation of the 2-D finite-difference time-domain code using high level shader language,” Applied Computational Electromagnetics Society Journal, vol. 23, no. 4, pp. 309–316, 2008.

A. Dziekonski, A. Lamecki, and M. Mrozowski, “GPU acceleration of multilevel solvers for analysis of microwave components with finite element method,” IEEE Microwave and Wireless Components Letters, vol. 21, no. 1, pp. 1–3, 2011.

M. Livesey, J. F. Stack, F. Costen, T. Nanri, N. Nakashima, and S. Fujino, “Development of a CUDA implementation of the 3D FDTD method,” IEEE Antennas and Propagation Magazine, vol. 54, no. 5, pp. 186–195, 2012.

R. Shao, D. Linton, I. Spence, P. Milligan, and N. Zheng, “Detection and GPU acceleration of 3D FDTD algorithms based on memory access patterns,” in Proceedings of the 2013 International Conference on Mechatronic Sciences, Electric Engineering and Computer, MEC 2013, pp. 2520–2526, China, December 2013.

T. T. Zygiridis, N. V. Kantartzis, and T. D. Tsiboukis, “GPU-accelerated efficient implementation of FDTD methods with optimum time-step selection,” IEEE Transactions on Magnetics, vol. 50, no. 2, 2014.

J. Francés, B. Otero, S. Bleda et al., “Multi-GPU and multi-CPU accelerated FDTD scheme for vibroacoustic applications,” Computer Physics Communications, vol. 191, pp. 43–51, 2015.

S. Nanthakumar, T. Lahmer, X. Zhuang, G. Zi, and T. Rabczuk, “Detection of material interfaces using a regularized level set method in piezoelectric structures,” Inverse Problems in Science and Engineering, pp. 1–24, 2015.

M. Cicuttin, L. Codecasa, B. Kapidani, R. Specogna, and F. Trevisan, “GPU accelerated time-domain discrete geometric approach method for maxwell's equations on tetrahedral grids,” IEEE Transactions on Magnetics, vol. 54, no. 3, 2018.

H. Ghasemi, H. S. Park, and T. Rabczuk, “A multi-material level set-based topology optimization of flexoelectric composites,” Computer Methods Applied Mechanics and Engineering, vol. 332, pp. 47–62, 2018.

T. M. Millington and N. J. Cassidy, “Optimising GPR modelling: a practical, multi-threaded approach to 3D FDTD numerical modelling,” Journal of Computers & Geosciences, vol. 36, no. 9, pp. 1135–1144, 2010.

I. Giannakis, A. Giannopoulos, and C. Warren, “Realistic FDTD GPR antenna models optimized using a novel linear/nonlinear full-waveform inversion,” IEEE Transactions on Geoscience and Remote Sensing, vol. 57, no. 3, pp. 1768–1778, 2019.

R. Luebbers, D. Ryan, and J. Beggs, “A two-dimensional time-domain near-zone to far-zone transformation,” IEEE Transactions on Antennas and Propagation, vol. 40, no. 7, pp. 848–851.

G. Mur, “Absorbing boundary conditions for the finite-difference approximation of the time-domain electromagnetic-field equations,” IEEE Transactions on Electromagnetic Compatibility, vol. 23, no. 4, pp. 377–382, 1981.

X. Wang, S. Liu, X. Li, and S. Zhong, “GPU-accelerated finite-difference time-domain method for dielectric media based on CUDA,” International Journal of RF and Microwave Computer-Aided Engineering, vol. 26, no. 6, pp. 512–518, 2016.