Mathematical Problems in Engineering

Volume 2014 (2014), Article ID 942106, 9 pages

http://dx.doi.org/10.1155/2014/942106

Research Article

Electromagnetic Computation and Visualization of Transmission Particle Model and Its Simulation Based on GPU

Yingnian Wu,1,2 Lin Zhang,¹ and Lan Mu¹

¹School of Automation Science and Electrical Engineering, Beihang University, Beijing 100191, China

²Beijing Key Laboratory of High Dynamic Navigation Technology, University of Beijing Information Science & Technology, Beijing 100192, China

Received 9 March 2014; Revised 4 July 2014; Accepted 18 July 2014; Published 7 August 2014

Academic Editor: Efstratios Tzirtzilakis

Copyright © 2014 Yingnian Wu 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.

Electromagnetic calculation plays an important role in both military and civic fields. Some methods and models proposed for calculation of electromagnetic wave propagation in a large range bring heavy burden in CPU computation and also require huge amount of memory. Using the GPU to accelerate computation and visualization can reduce the computational burden on the CPU. Based on forward ray-tracing method, a transmission particle model (TPM) for calculating electromagnetic field is presented to combine the particle method. The movement of a particle obeys the principle of the propagation of electromagnetic wave, and then the particle distribution density in space reflects the electromagnetic distribution status. The algorithm with particle transmission, movement, reflection, and diffraction is described in detail. Since the particles in TPM are completely independent, it is very suitable for the parallel computing based on GPU. Deduction verification of TPM with the electric dipole antenna as the transmission source is conducted to prove that the particle movement itself represents the variation of electromagnetic field intensity caused by diffusion. Finally, the simulation comparisons are made against the forward and backward ray-tracing methods. The simulation results verified the effectiveness of the proposed method.

1. Introduction

Electromagnetic calculation plays an important role in both military and civic fields. It is realized that simulation of an electromagnetic environment is useful and can give much information in communication systems design and radar testing. It is of great significance to calculate the distribution status of the electromagnetic field.

Some methods and models are proposed for calculation of indoor and outdoor wave propagation, that is, methods based on the numerical solution of Maxwell’s equations such as finite-difference time-domain method (FDTD) [1] and parabolic equation method (PE) [2–6], methods based on geometrical optics principle such as ray-tracing (RT) [7–10] method, and so forth.

To calculate electromagnetic wave in a large range, FDTD method brings heavy burden in calculation and also requires huge amount of memory [11–14]. PE method performs well for complex terrain but has significant limitation in the calculation of propagation angle. Generally speaking, RT method is applicable to the electromagnetic calculation in a large range.

Ray-tracing method is based on geometrical optics (GO) principle, which decides reflection or refraction through simulating the propagation route of rays [15–17]. It also introduces diffraction ray in the case of obstacles. According to GO theory, improved geometry theory of diffraction (GTD) [18] and uniform theory of diffraction (UTD) [19–21] are established. Ray-tracking is classified into forward ray-tracking and backward ray-tracking [22, 23]. Backward ray-tracking is generally used because forward ray-tracking has larger error. But backward ray-tracking has huge calculation when calculating space field, because every point in the space has to be recalculated.

Particle method is used to simulate the macrostate based on a huge amount of individual particle which has a process of generation, motion, and vanishing.

Based on forward ray-tracing method and particle method, a Transmission Particle Model for calculating electromagnetic distribution status in space is presented in this paper. Using the GPU to accelerate computation and visualization can reduce the computational burden on the CPU. The movement of a particle obeys the principle of the propagation of electromagnetic wave, and then the particle distribution density in space reflects the distribution status of electromagnetic distribution status. Since the particles in TPM are completely independent, it is very suitable for the parallel computing based on GPU. Updated GPU particle data is stored in display memory, so the data does not need to be transferred back to the memory when making visualization of the particle, and it can directly draw the particle movement after transmission.

In this paper, we present Transmission Particle Model and use GPU based on CUDA, realize particles transmit, updates, demise, and other operations, and achieve particle model computation and synchronization visualization.

2. Description of TPM

The basic idea of TPM is that particles are generated and transmitted intermittently in the origin of electromagnetic field, and the electromagnetic distribution status in space is obtained through calculating the movement of each particle.

So far there are many methods for electromagnetic calculation, such as finite-difference time-domain method (FDTD), parabolic equation method (PE), and ray-tracing method (RT). To calculate electromagnetic wave in a large range, FDTD method brings heavy burden in calculation and also requires huge amount of memory. PE method performs well for complex terrain but has significant limitation in the calculation of propagation angle. Generally speaking, RT method is applicable to the electromagnetic calculation in a large range.

2.1. Particle Transmission

From the generation in the origin of electromagnetic field, the particle transmission density mainly depends on antenna’s angle factor . and are the angle components in sphere coordinates.

Take the origin of electromagnetic field as a vertex; partition the space with equal interval in all directions. The sampling intervals for and are and , respectively; usually we have .

Then we obtain

Denote as the particle numbers in one transmission in the region , , whose maximum value is ( is the product of and the product of and , , and are given constants; if , then we have , ):

The initial coordinates are , initial speed is , initial distance , the initial interval angle with other adjacent particles is , the initial life time is , and the initial particle status is active, where

is the distance between the initial position of particle transmission and the electromagnetic transmission source, which is supposed to be the origin of coordinate.

2.2. Particle Movement

The movement of particle is dynamic; that is, with the elapse of time, the current properties, such as position coordinates, speed , movement distance , angle , life time , and state (active or dead), will affect the properties of the particles in the particle list.

In Figure 1, is the total movement distance of the particle, is the distance between particles in the same batch, and is the distance particle of two consecutive batches within time interval . Particles are consecutively transmitted with equal time interval .

Figure 1: The distance between the particles.

Since is very small, at initial time, we have

At time , we have

Since is very small, we have . Consequently, particles are approximately in uniform distribution in space; particle density in can represent the energy of the electromagnetic energy. In the next section we will give the expression between particle density and Poynting vector mean value for electric dipole.

The movement speed of particle is

Correspondingly, we have

2.3. Particle Reflection

When electromagnetic particles come to the interface of two media, they reflect according to reflection theory, some of them may vanish owing to reflection, and the probability of vanishing is decided by reflection factor .

The liveability after reflection is

Then the vanishing probability is

If the electric field intensity vector is parallel to incidence plane (parallel polarization), the reflection factor of uniform plane wave in [21] is

If the electric field intensity vector is vertical to incidence plane (vertical polarization), the reflection factor of uniform plane wave in [10] is

After reflection, some particles vanish. The transmission interval angle between survived particles and their adjacent particles has been changed to :

2.4. Particle Diffraction

In Geometry Optics, there are three kinds of diffractions, that is, edge diffraction, spire diffraction, and surface diffraction. Since the attenuation of spire diffraction and surface diffraction is very fast, we can ignore the two kinds of diffractions and only take edge diffraction into account.

When the particle diffraction occurs, the position of the particle needs to be modified according to the diffraction principle of electromagnetic wave. And further, the death probability of the particle is calculated based on the diffraction factor, which can be used to decide the death of the particle.

Given a distance , there is particle diffraction within the range of around the edge of terrain. After the diffraction, the particle is located on the surface of a cone whose vertex is the diffraction point, the axis is the tangent line of the contour around the diffraction point, and the angle between the generatrix and the axis is equal to the angle between the incidence particle and the tangent line.

The liveability after diffraction is

Then the vanishing probability after diffraction is

Similar to reflection, the transmission interval angle between the survived particles after diffraction, that is, , will be changed to :

3. Model Verification

Electric dipole is the most simple electromagnetic transmission antenna. This paper is focused on the analysis of the electromagnetic field generated by electric dipole.

3.1. Electric Dipole

In unbounded ideal media, electromagnetic wave propagates without reflection, refraction, scattering, and absorption; there is only energy loss caused by diffusion. In calculation of electromagnetic transmission in a large area, the average power density of close field is near zero, that is, no electromagnetic power outputs in the close field. Thus we only need to do calculations for the far-reaching field .

In the far-reaching electromagnetic field, the electric field vector and magnetic field vector can be expressed approximately in [10] as follows:

is the virtual value of current, is the length of electric dipole, is the number of waves, ( is the wavelength of electromagnetic wave), is the angular frequency, and is the dielectric constant in free space.

Poynting vector mean value of any point in the space is

The radiation field generated by electric dipole is not only relevant to the distance from field point to the source but also relevant to the and in spherical coordinates.

By (24) and (25) can be seen, electric field strength and magnetic field strength of electric dipole changes with the angle changing, which can be expressed as a function called antenna’s angle factor.

The antenna’s angle factor of electric dipole is

According to function , we can get the angle factor figure of electric dipole which is shown in Figure 2.

Figure 2: Angle factor figure of electric dipole.

3.2. Verification Deduction of Electromagnetic TPM

Consider the coordinate origin is the transmission source with sampling intervals and , ; then after long enough time, the particle density at a point is within the range: , . Then the number of transmitted particles at the centre of the above-mentioned area is , where

The distance range of these particles to the transmission source is :

The distance between the transmission source and the particles transmitted at the previous time instant is

where

Since these particles are uniformly distributed in local area, then we have , , .

The volume of the area is

where

Substituting (27) into (26), we get

Then the density of particles is

Considering (21), (22), and (25), then we obtain

Since , then , and .

Further we have

On account of , , , then

In order to express Poynting vector mean value by the density of particles and since the antenna’s angle factor of electric dipole is

we have

Then we can get

Substituting this into (17), we obtain

where

Consider is a constant, and then we know that the electromagnetic status in uniform space can be calculated through TPM. Poynting vector mean value of any point in the space is identical to the particle density of that point.

4. Computation and Synchronization Visualization Based on GPU

4.1. GPU Acceleration Technology

In recent years, due to the limit of the integrated circuit components, we encountered bottlenecks when we want to enhance the central processing unit (CPU) clock frequency on existing infrastructure. When upgrading the performance of a single processor core becomes difficult, multicore CPU was launched, and people are turning to the Graphics Processing Units (GPU). In the past GPU was primarily used for graphics rendering and display, but with the improvements of GPU performance and functional structure, GPU is not only used for graphics processing and rendering, but also starts to become a powerful assistant for CPU computing.

The current mainstream GPU has the following characteristic features: (1) supporting IEEE 32-bit floating-point arithmetic and output; (2) providing flexible programmability, dynamic branching, and loop control functions; (3) supporting multiple rendering, avoiding multiple data transfer between the CPU and GPU; and (4) supporting texture access function. It can store data to the texture easily and get index access. (5) It supports output data to the texture features, which can improve data output speed and reduce the overhead of data readback. Using the GPU to accelerate computation and rendering can reduce the computational burden on the CPU, free CPU from the heavy computing tasks, and meet the requirements of other computing.

As the inventor of the GPU, NVIDIA launched a ​​common computing programming environment CUDA, completely shielding the underlying graphics API, which makes GPU programming almost transparent for developers and lays the foundation for wide use of GPU common computing.

4.2. TPM Computation Based on CUDA

The process for the calculation of electromagnetic distribution status based on CUDA with the particle transmission model can be described as follows (Figure 3).

Figure 3: The process for the calculation of electromagnetic based on CUDA.

Step 1. Transmit particles, initiate each of the particles, and put the transmitted particles to the particle list.

Step 2. Update the properties of the particles in the list for the next time instant.

Step 3. Judge if the particle is below the ground. If so, reflection or diffraction is needed. Next, the choice between reflection and diffraction is made based on whether the particle is on the edge of spire of the terrain.

Step 4. If reflection is needed, modify the position of the particle according to the reflection principle of electromagnetic wave. And further, calculate the death probability of the particle based on the reflection factor and decide the death of the particle based on the vanishing probability. If diffraction is needed, the process is similar to that of particle reflection.

Step 5. Judge if the particle is outside the boundary, and then eliminate the outside particles and the death particles from the list.

Step 6. If the number of the listed particles is not stable (i.e., the particles within the calculation area have not reached equilibrium), then return to Steps 1–5.

Step 7. If the number of the particles in the list is stable, partition the space into grids and then count and output the number of particles within each grid.

4.3. Computation and Visualization of Particle Model of Synchronization

Since the particle in Transmission Particle Model is completely independent, TPM is very suitable for the parallel computing based on CUDA. Updated GPU particle data is stored in display memory, so the data does not need to be transferred back to the memory when making visualization of the particle. It can directly draw the particle movement after transmission.

In order to achieve the results, we use CUDA to compute and draw directly on the display memory. As we need to transfer data between OpenGL and CUDA, we should establish buffer zones between these two groups of API for the position coordinates and speed. We can get the actual address of the device memory through the kernel function. Firstly, call cudaGraphicsMapResources() function to tell CUDA runtime mapping the OpenGL shared memory, and then call cudaGraphicsResourceGetMappedPointer() to request a pointer to the mapped resource. Thus, the OpenGL can share with CUDA buffer.

5. Simulation and Analysis

5.1. Comparison of TPM with RT

Electromagnetic particle model is based on ray-tracing method, but it is rather different from ray-tracing. In ray-tracing method, the space distances between transmission rays are equal; we calculate the initial density on different directions and need to calculate the variation of the density of electromagnetic field caused by diffusion in space. As to electromagnetic particle model, the ray is replaced by particles, the energy of each particle is a constant, the intensity of electromagnetic field is represented by the density of particles, the transmission density on different directions depends on the angle factor of the antenna, and the movement of particle accords with the attenuation of diffusion in space. Consequently, the particle movement itself represents the variation of electromagnetic field intensity caused by diffusion.

To calculate the electromagnetic field intensity at a point, forward ray-tracing algorithm adopts the method of receiving sphere to define the receiving radius and count the numbers of rays in the scope of the receiving sphere. But the definition of receiving sphere may bring in significant errors. Consider there is only direct transmission in an infinite space; one receiving sphere should receive only one ray. The receiving radius equals the distance between two rays around the point. If the receiving sphere is defined larger than the distance, there could be two rays in the receiving sphere; if the receiving sphere is defined smaller than the distance, there could be no rays in the receiving sphere. That is the reason for errors in the definition of the receiving sphere. The energy of each particle in electromagnetic particle model is less than that of a ray in ray-tracing method, thus yielding the smaller error. Of course, in this way the deficiency of significant prediction errors in forward ray-tracing method can be overcome.

Compared with backward ray-tracing, TPM can calculate the electromagnetic distribution status in the whole area at one time, but the backward ray-tracing needs to calculate each point in the area once again. In other words, TPM needs less calculation than backward ray-tracing to get electromagnetic distribution status.

To verify the effectiveness of the proposed TPM, we made a comparison between the proposed method and the ray-tracing method with a simplified terrain, and the transmission source is an electric dipole antenna in the middle of the area.

In order to compare TPM with ray-tracing method, we choose several points to calculate the propagation loss according to the ray-tracing method and the particle method, respectively. The height for calculation is 100 m. Because backward ray-tracing is commonly used, the result is calculated by backward ray-tracing as standard. The results are shown in Figure 4. The beginning of the curve is not drawn because this part is under the ground. There are errors in the calculations using the forward ray-tracing method. We can see that the power is bigger than standard near 95 m because the receiving sphere is defined larger in Figure 4(a), and the receiving sphere is defined lower near 78 m in Figure 4(b). The calculation results of TPM and backward ray-tracing are almost the same; thus TPM is effective for electromagnetic calculation and can overcome the deficiency of errors in forward ray-tracing method.

(a)

(b)

Figure 4: Compare TPM with ray-tracing method.

5.2. Comparison of Transmission Particle Model and Wireless InSite

Wireless InSite is complex electromagnetic environment predictive analysis simulation commercial software. The software is based on UTD/GTD theory and mainly uses ray-tracing method combined with the improved FDTD method to calculate the electromagnetic propagation. The calculation model includes urban canyon model, fast 3D city models, full three-dimensional model, the vertical plane model, the free-space model, urban canyon FDTD model, FDTD model with sliding windows, HATA model, COST-HATA model, and WI real-time models. By combining electricity field with the antenna pattern, it calculates the path loss, time of arrival, angle of arrival of the electromagnetic field, and the electromagnetic trend distribution.

To further verify the feasibility of the model of Transmission Particle Model, we calculated the electromagnetic distribution in the same terrain using particle transmission simulation model and Wireless InSite and displayed the electromagnetic attenuation trend extracting surfaces 5 meter far from ground. Figure 5(a) shows the calculation results of Wireless InSite, which set 20 20 receivers, and each colour block represents the energy for one receiver. Figure 5(b) shows the calculation of Transmission Particle Model. The space is divided into 100 100 100 grids. Figure 5(b) is the extracted 100 100 grids 5 meters from the ground constituting a surface. As we can see from the comparison of the two figures, the electromagnetic attenuation of Transmission Particle Model and Wireless InSite is basically in the same situation, which verified the reliability of the Transmission Particle Model.

(a) Wireless InSite

(b) Transmission Particle Model

Figure 5: Visualization comparing TPM with Wireless InSite software (dBm).

The calculation time of TPM is less than the WI using RT method. The calculation time comparison is shown in Figure 6.

Figure 6: Calculation time comparing TPM with Wireless InSite (RT).

When Wireless InSite software needs to calculate the electromagnetic field strength of a point, it uses the method of setting the receiver and calculates the received electromagnetic energy for the receiver. If we want to calculate the electromagnetic field distribution trend of the entire space, we need to provide intensive receiver in the space at equal intervals. Transmission Particle Model can calculate electromagnetic distribution trend in the entire space of the region. In the calculation of the electromagnetic distribution trend in the entire space, the calculation amount of Transmission Particle Model is less than that of Wireless InSite software.

6. Conclusion

Based on ray-tracing method and the particle method, this paper constructs TPM for electromagnetic field calculation. The details, including the initiation of the particles, updating the properties of the particles, the reflection and diffraction calculation, are also presented. Since the particles in TPM are completely independent, we used parallel computing based on GPU for the computation. Updated GPU particle data is stored in display memory, when making visualization of the particle, and it can directly draw the particle movement after transmission. Deduction verification of TPM is conducted based on the electromagnetic field generated by electric dipole. Finally, the simulation comparisons are made against the ray-tracing method and the commercial software WI. The simulation results verified the effectiveness of the proposed method. In summary, TPM is an effective tool for electromagnetic calculation, such as the electromagnetic distribution status in the space with complex terrain.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This work is supported by Scientific Research Common Program of Beijing Municipal Commission of Education (KM201411232007), Open Research Project of the Beijing Key Laboratory of High Dynamic Navigation Technology (Grant no. HDN2014007), and the National Natural Science Foundation of China (Grant no. 61031001), Beijing Natural Science Foundation (4142017).

References

J. W. Schuster and R. J. Luebbers, “Comparison of GTD and FDTD predictions for UHF radio wave propagation in a simple outdoor urban environment,” in Proceedings of the IEEE Antennas and Propagation Society International Symposium, pp. 2022–2024, Digest, Montreal, Canada, July 1997.

C. Kuifu and Z. Senwe, “DFT accelaration and its application in non-stationary random vibration calculation,” Journal of Vibration Engineering, vol. 12, pp. 91–96, 1999 (Chinese).

M. F. Levy, Parabolic Equation Methods for Electromagnetic Wave Propagation, Institution of Engineering and Technology Press, London, UK, 2000.

Z. Senwen and C. Kuifu, “The monment equation method for multi-degree of freedom system with strongly nonlinear coupling and parametric excitation,” Acta Mechnica Sinica, vol. 25, pp. 362–368, 1993 (Chinese).

S. Amiri and S. M. Hosseini, “Second-order method for solving 2D nonlinear parabolic differential equations based on ADI method,” International Journal of Modeling, Simulation, and Scientific Computing, vol. 1, no. 1, pp. 133–146, 2010.

V. P. Tanana, “An order-optimal method of solving an inverse problem for a parabolic equation,” Journal of Applied and Industrial Mathematics, vol. 3, pp. 395–400, 2009.

Y. Cocheril and R. Vauzelle, “A new ray-tracing based wave propagation model including rough surfaces scattering,” Progress in Electromagnetics Research, vol. 75, pp. 357–381, 2007.

B. R. Epstein and D. L. Rhodes, “GPU-accelerated ray tracing for electromagnetic propagation analysis,” in Proceedings of the IEEE International Conference on Wireless Information Technology and Systems, (ICWITS '10), pp. 1–4, September 2010.

P. G. Brown and C. C. Constantinou, “Investigations on the prediction of radiowave propagation in urban microcell environments using ray-,” IEE Proceedings: Microwaves, Antennas and Propagation, vol. 143, no. 1, p. 36, 1996.

T. E. Athanaileas, G. E. Athanasiadou, G. V. Tsoulos, and D. I. Kaklamani, “Parallel radio-wave propagation modeling with image-based ray tracing techniques,” Parallel Computing, vol. 36, no. 12, pp. 679–695, 2010.

S. Aono, M. Unno, and H. Asai, “Alternating-direction explicit FDTD method for three-dimensional full-wave simulation,” in Proceedings of the 60th Electronic Components and Technology Conference (ECTC '10), pp. 375–380, June 2010.

N. Endoh, T. Tsuchiya, and S. Matsumoto, “Analysis of acoustic characteristics of aero ultrasonic sensor calculated by 3D-FDTD,” in Forum Acusticum Budapest 2005: 4th European Congress on Acustics, Scientific Society for Optics, Budapest, Hungary, 2005.

J. E. Bendz, H. G. Fernandes, and M. K. Zuffo, “Fast and accurate finite-difference time-domain method for large three-dimensional simulations,” in Proceedings of the 19th IASTED International Conference on Modelling and Simulation, pp. 20–25, Quebec City, Canada, May 2008.

F. Yalcin Yamaner and A. Bozkurt, “P5E-4 3D transient analysis of ultrasound propagation using finite difference time domain method and its experimental verification,” in IEEE Ultrasonics Symposium, pp. 2295–2298, 2007.

G. E. Athanasiadou and A. R. Nix, “A novel 3-D indoor ray-tracing propagation model: the path generator and evaluation of narrow-band and wide-band predictions,” IEEE Transactions on Vehicular Technology, vol. 49, no. 4, pp. 1152–1168, 2000.

E. G. Papkelis, I. Psarros, I. C. Ouranos et al., “A radio-coverage prediction model in wireless communication systems based on physical optics and the physical theory of diffraction,” IEEE Antennas and Propagation Magazine, vol. 49, no. 2, pp. 156–165, 2007.

A. Escobar, L. Lozano, H. Cadavid, and M. F. Cátedra, “A new ray-tracing acceleration technique for radio propagation,” in Proceeding of the IEEE Antennas and Propagation Society International Symposium (APSURSI '10), pp. 1–4, Toronto, Canada, July 2010.

J. B. Keller, “Geometrical theory of diffraction,” Journal of the Optical Society of America, vol. 52, pp. 116–130, 1962.

R. G. Kouyoumjian and P. H. Pathak, “A uniform geometrical theory of diffraction for an edge in a perfectly conducting surface,” Proceedings of the IEEE, vol. 62, no. 11, pp. 1448–1461, 1974.

M. B. Tabakcioglu and A. Kara, “Comparison of improved slope uniform theory of diffraction with some geometrical optic and physical optic methods for multiple building diffractions,” Electromagnetics, vol. 29, no. 4, pp. 303–320, 2009.

G. Z. Ni, Engineering Principle of Electromagnetic Field, Higher Education Press, Beijing, China, 2009.

Y. Corre and Y. Lostanlen, “Three-dimensional urban EM wave propagation model for radio network planning and optimization over large areas,” IEEE Transactions on Vehicular Technology, vol. 58, no. 7, pp. 3112–3123, 2009.

S. Priebe, M. Jacob, C. Jastrow, T. Kleine-Ostmann, T. Schrader, and T. Kürner, “A comparison of indoor channel measurements and ray tracing simulations at 300 GHz,” in Proceedings of the 35th International Conference on Infrared, Millimeter, and Terahertz Waves (IRMMW-THz '10), pp. 1–2, September 2010.