# Analysis of Electromagnetic Propagation from MHz to THz with a Memory-Optimised CPML-FDTD Algorithm

Authors: A. Rodríguez-Sánchez; C. Couder-Castañeda; J. J. Hernández-Gómez; I. Medina; S. Peña-Ruiz; J. Sosa-Pedroza; M. A. Enciso-Aguilar
Journal: International Journal of Antennas and Propagation
 (2018)

---

## Abstract

FDTD method opened a fertile research area on the numerical analysis of electromagnetic phenomena under a wide range of media and propagation conditions, providing an extensive analysis of electromagnetic behaviour like propagation, reflection, refraction, and multitrajectory phenomena. In this paper, we present an optimised FDTD-CPML algorithm, focused in saving memory while increasing the performance of the algorithm. We particularly implement FDTD-CPML method at high frequency bands, used in several telecommunications applications as well as in nanoelectromagnetism. We show an analysis of the performance of the algorithm in single and double precision, as well as a stability of the algorithm analysis, from where we conclude that the implemented CPML ABC constitutes a robust choice in terms of precision and accuracy for the high frequencies herein considered. It is important to recall that the CPML ABC parameters provided in this paper are fixed for the tested range of frequencies, that is, from MHz to THz.


---

## Body

## 1. Introduction

In the last 30 years, a wide range of absorption techniques in computational boundaries has been developed for FDTD as well as for other methods as finite and spectral elements. The main purpose of absorbing boundary conditions (ABCs) is to avoid spurious reflection of energy within the physical domain, at any incidence angle and operation frequencies, considering the thinner possible absorption zone. Boundary conditions for electromagnetic phenomena simulations are absolutely important as they allow to truncate the physical domain in the regions of interest, to simulate an infinite region without spending infinite computational resources (as RAM memory). A good choice of ABCs involves three key features: easiness of implementation, efficiency of the algorithm, and precision for long time periods. Actually, mostsolvers in the FDTD field implement PML ABCs [1], due to its extremely implementation easiness. Nevertheless, Berenger’s method heavily relies on the splitting of the fields in the boundary region, so for large domains, it constitutes a computationally expensive method. Convolutional PML [2] alleviates this feature due to the fact that it does not split any field in the boundary region, so the number of memory arrays for it to be implemented is reduced with respect to Berenger’s formulation.Besides, Berenger’s PML ABC [1] is not fully capable of avoiding all spurious reflections [3, 4]. CPML ABC introduced by Komatitsch and Martin [3] is a convenient version of the convolutional PML for its implementation in different computational architectures, because it reduces the RAM memory usage in absorption regions, while providing an excellent precision. More recently, Martin and Couder-Castañeda introduced an unsplit convolutional perfectly matched layer (CPML) technique to efficiently absorb the compressible viscous flow without causing any numerical instability over long periods of time [5].Nevertheless, most of the research done in the numerical propagation of electromagnetic phenomena with FDTD-CPML method has been focused on providing CPML ABC parameters to effectively absorb electromagnetic energy for relatively low frequencies (from Hz to MHz). Actually, FDTD studies for high and ultrahigh frequencies pose numerical difficulties as the memory and processing times required to resolve accurately the problem for such frequencies are still much higher than those provided by current computers, particularly when propagating in dispersive media [6, 7]. This problem has derived in two research directions. The first one has tried to use the recent progress in high performance parallel computing, such as those provided by the development of GPUs, [8–11] as well as multiprocessing coprocessors like Xeon Phi [12–16]. The other branch constitutes the extreme optimisation of current FDTD-CPML algorithms to be able to be successfully executed in current computers. This work belongs to the latter research branch.In this paper, we particularly tackle the problem of drastically reducing the memory requirements at the absorption boundaries by optimising the CPML ABC developed by Komatitsch and Martin [3] for the Maxwell equations in FDTD method, in order to study its precision and accuracy in different propagation scenarios for the electric transverse propagation mode. We determine a set of CPML ABC parameters which effectively absorb all the energy in the boundary region for the whole frequency range from MHz to THz, which as to our knowledge has not been previously reported in literature. We tested the FDTD-CPML implemented method in single and double precision, in two different implementations of the algorithm, in order to minimise the memory usage as well as computing time.This paper is organised as follows: in Section2, we present the electromagnetic theory that gives rise to the FDTD method herein implemented. In Section 3, we introduce the optimised CPML ABC formulation we use, while in Section 4, we present numerical experiments: electromagnetic propagation in free space, with a parabolic reflector, in a coplanar nanowaveguide, as well as in common office indoors, along with their respective error analysis. In Section 5, we provide a performance study of the memory usage with double and single precision, as well as of computing time. Finally, in Section 6, we condense the results herein obtained.


## 2. Propagation Equations and FDTD Method

Electromagnetic phenomena are governed by Maxwell equations, which for a linear, homogeneous, and isotropic media, in the international system of units, are given by [17]
(1)∇→⋅E→=1ερs,(2)∇→⋅H→=0,(3)∇→×E→+μ∂H→∂t=0,(4)∇→×H→−∈∂E→∂t=J→s+σE→,where E→ and H→ are the electric and magnetic fields, respectively, ε and μ are permittivity and permeability constants that depend on the medium in which the electromagnetic field is propagated, and σ is the electric conductivity of the material. It is necessary to remark that the speed of light in such a medium is given by c=1/εμ. The external current and electric charge densities, J→s and ρs, respectively, are the sources of electromagnetic field, and they can be functions of space and time, that is, J→s=J→sr→,t and ρs=ρsr→,t.For the sake of simplicity, the system of (1), (2), (3), and (4) can be reduced to a bidimensional system, for instance to an electromagnetic propagation in x,y plane, by confining the energy to propagate in such a plane (i.e., the Poynting vector). The resulting equations for the transverse electric mode (TE) in the absence of sources, and considering Cartesian coordinates, are thus
(5)∂Ez∂t=1ε∂Hy∂x−∂Hx∂y−σEz,(6)∂Hy∂t=1μ∂Ez∂x,(7)∂Hx∂t=−1μ∂Ez∂y.In order to cast (5), (6), and (7) in a second order in time and space FDTD method, the index notation for the considered staggered mesh shown in Figure 1 is followed. If time evolution is discretised with the index n, then the finite differences version of such equations are, respectively,
(8)Ezi,jn+1=Ezi,jn1−Δtσi,j/2εi,j1+Δtσzi,j/2εi,j+Δtεi,j1+Δtσi,j/2εi,jHyi,j+12n+12−Hyi,j−12n+12Δx−Δtεi,j1+Δtσi,j/2εi,jHxi+1/2,jn+1/2−Hxi−1/2,jn+1/2Δy,(9)Hyi,j+1/2n+1/2=Hyi,j+1/2n−1/21−Δtσi,j+1/2/2μi,j+1/21+Δtσi,j+1/2/2μi,j+1/2+Δtμ1+σi,j+1/2Δt/2μi,j+1/2Ezi,j+1n−Ezi,jnΔx,(10)Hxi+1/2,jn+1/2=Hxi+1/2,jn−1/21−Δtσi+1/2,j/2μi+1/2,j1+Δtσi+1/2,j/2μi+1/2,j+Δtμ1+σi+1/2,jΔt/2μi+1/2,jEzi+1,jn−Ezi−1,jnΔy.Figure 1
Indices notation for the considered mesh.


## 3. CPML ABC Formulation

A very suitable version of the CPML ABC is given in [18]; this is because its implementation in a FDTD code without a PML ABC is simply obtained by replacing the spatial derivative ∂x, with 1/k∂x+Ψ, where Ψ is the so-called convolutional term of the CPML ABC—whose time evolution is the same as the time evolution of the other variables—and k is the stretching function always greater than 1 [5] and given by (13).According to [18], Ψ can be implemented as a recursive sum in time, so we can update the memory variable of the F (H or E) field in the i=x,y direction for each time step n as
(11)ΨinF=biΨin−1F+ai∂Fn−1.For instance, the auxiliary equation for∂xHy takes the form
(12)Ψ∂xHyi,jn+1/2=bxiΨ∂xHyi,jn+1/2+axi∂xHyi,jn+1/2.The vectorskx, bx, and ax are estimated, within the absorption region, as
(13)kxq=1+kmax−1.0xqnormm,where xqnorm = Th−qΔx/Th, Th is the CPML thickness, q=0,1,2…N, and N is the number of layers of the CPML ABC. It is important to note that within the physical domain, bx=0, cx=0, and kx=1; within the CPML
(14)bxq=e−dxq/kxq+axqΔt/ε0,(15)axq=dxqbxq−1dxqkxq+kxq2axq,where ε0 is the vacuum permittivity,
(16)dxq=dmaxxqnormmis the damping function profile, dmax=0.750.8m+1/Δxμ0/ε0,
(17)axq=amax1.0−xqnormmais the shift in the frequency domain [5], and m and ma are the degree of polynomial variation for dx and ax, respectively.The vectorskxh, bxh, and axh are calculated in the intermediate steps j+12, while the vectors kyh, byh, and ayh are calculated in the intermediate steps i+1/2.Equations (8), (9), and (10) with the implemented CPML take the form
(18)Ezi,jn+1=Ezi,jn1−Δtσi,j/2εi,j1+Δtσzi,j/2εi,j+Δtεi,j1+Δtσi,j/2εi,j1kxiHyi,j+1/2n+1/2−Hyi,j−1/2n+1/2Δx+Ψ∂xHyi,jn+1/2−Δtεi,j1+Δtσi,j/2εi,j1kyjHxi+1/2,jn+1/2−Hxi−1/2,jn+1/2Δy+Ψ∂yHxi,jn+1/2,(19)Hyi,j+1/2n+1/2=Hyi,j+1/2n−1/21−Δtσi,j+1/2/2μi,j+1/21+Δtσi,j+1/2/2μi,j+1/2+Δtμ1+σi,j+1/2Δt/2μi,j+1/21kxhiEzi,j+1n−Ezi,jnΔx+Ψ∂xEzi,j+1/2n,(20)Hxi+1/2,jn+1/2=Hxi+1/2,jn−1/21−Δtσi+1/2,j/2μi+1/2,j1+Δtσi+1/2,j/2μi+1/2,j+Δtμ1+σi+1/2,jΔt/2μi+1/2,j1kyhjEzi+1,jn−Ezi−1,jnΔy+Ψ∂yEzi+1/2,jn.Equations (18), (19), and (20) govern the applications herein treated.


## 4. Numerical Experiments

We considered four different propagation scenarios:(i)
Free space propagation
(ii)
Parabolic reflector
(iii)
Coplanar nanowaveguides (CNWG)
(iv)
Indoor propagationTo corroborate the quality of the numerical solution, it is necessary to determine the reflection error due to the presence of the CPML ABC, so reference solutions must be obtained for each of the applications previously mentioned. These reference solutions consist in performing the respective electromagnetic propagation with an extended physical domain, not allowing the waves to reach the boundaries of the extended domain to avoid spurious reflected waves to propagate within the extended domain.An advantage of the usage of the presented CPML ABC consists in a drastic reduction of the number of memory arrays in the bidimensional algorithm [4], so it can be easily implemented for GPU processing, as those cards possess a limited amount of memory. The CPML ABC for the FDTD method requires to allocate the value of the time derivatives in a memory variable, which is implemented in two ways in this paper: the first, by allocating the memory variables for all the domain, and the second, by allocating them only in the absorption region (see Figure 2). Clearly, by allocating them just in the CPML region, we save memory, but with an increasing complexity in its implementation. By allocating the memory variables in all the domain, we require three nested cycles for calculation, but those cycles might be divided to save memory, producing eleven additional cycles to reach a total of fourteen nested cycles, which might decrease the performance of the vectorisation processes [19, 20].Figure 2
CPML ABC implementation options.The novel set of CPML parameters used for the herein reported numerical experiments are shown in Table1. These parameters have shown to be optimised for the whole range of frequencies tested in this paper. They were obtained by successive numerical experiments as follows: clearly, the convolutional term Ψ of the CPML given by (11) fully depends on the vectors ax and bx ((14) and (15)), which in turn fully depends on kx, dx and ax ((13), (16), and (17), respectively,). Thus, the parameters to be adjusted to study the appropriateness and efficiency of the CPML ABC are those involved in calculating kx, dx, and ax, which from (13), (16), and (17), turn to be kmax, m, amax, ma, and N. In order to find the values of the parameters reported in Table 1, it is better to start from those reported in [2]. With those parameters fixed, we first varied the parameters regarding the vector kx for the range of frequencies herein considered until the parameters which do not distort field wavefronts are found. Then, we varied the parameters regarding the vector ax for all the range of frequencies until those which better absorb the incident fields are found.Table 1
Parameter values of the implemented CPML ABC.
Parameter
Value
μ0
4π10−7
ε0
1c2μ0
m
3.0
ma
1.0
kmax
3.0
amax
0.08
N
20The FDTD-CPML method was implemented in FORTRAN; beside its high speed for numerical algorithms, it eases the implementation because indexes might take positive or negative values, so they can start in (0,0) in the physical domain, so if the number of discrete points in thex direction is IMAX and in the y direction is JMAX, then the whole domain including the CPML might be labelled in the x direction with i=Nxs…Nxe and in the y direction with j=Nys…Nye, where Nxs=−N+1, Nxe=IMAX+N, Nys=−N+1, Nye=JMAX+N, and N is the number of layers of the CPML ABC. The index notation can be observed in Figure 1.Following the first implementation scheme, the flux diagram of the algorithm is shown in Figure3, while the respective flux diagram by allocating memory variables just in the absorption region is shown in Figure 4. We compare the performance in both implementations for each numerical experiment presented in this paper.Figure 3
Flux diagram of the numerical algorithm using memory variables in all the domain. The number of used cycles is numbered asV1, V2, V3, and V4.Figure 4
Flux diagram of the numerical algorithm using memory variables only in the absorption region. The cycles are numbered asM1⋯M14.The spatial incrementsΔx and Δy are estimated considering λ/20 (λ being the propagated wavelength), while Δt is estimated with the Courant stability criterion Δt=1/c1/Δx2+1/Δy2. The sources of electromagnetic field are modelled as sin2π×f×n×Δt, where f is the frequency and n is the time step number.For free space and parabolic reflector propagations, three operating frequencies were explored: 2.45, 5.00, and 20.00 GHz. The CNWG was tested at 100 THz while the indoor propagation was explored at a 2.45 GHz Wi-Fi frequency.For each numerical experiment, results are presented as follows: snapshots of field propagation are shown at different time steps for both the CPML and the reference solutions; for reasons of the paper’s length, we only show snapshots for an operating frequency of 2.45 GHz (except for the CNWG, which is of 100 THz), provided their behaviour for the other operating frequencies is similar. We also perform an analysis of the energy decay within the physical domain when the source is turned off (at the middle of the simulation time). The global energyU within the physical domain, at time step n, is given by
(21)Un=∑i=1IMAX∑j=1JMAXHxi,jnEzi,jn+Eyi,jnHzi,jn.To quantitatively verify the CPML performance, we set numerical viewers at the same positions in both the CPML and reference solutions for all the numerical experiments, to record the electromagnetic field in such points as a function of time. We then compare the recorded values of the viewers for both the reference and the ABC cases. For the reasons of the paper’s length as well as due to the likeness between them, only the viewer plots for 2.45 GHz are shown (100 THz for the CNWG).


### 4.1. Free Space Propagation

The physical domain set up is 1.25 m × 2.25 m while the reference (extended) domains are of 20 m × 20 m, 10 m × 10 m, and 5 m × 5 m for 20.00 GHz, 5.00 GHz, and 2.45 GHz, respectively, with the source located in the centre. Numerical viewer locations are shown in Figure5. The results of the simulation for the CPML and reference solutions are shown in Figure 6, where the absorption in the CPML region is observed.Figure 5
Configuration for the free space propagation scenario. Orange lines indicate the beginning of the absorption region (CPML); the orange dot in the centre indicates the location of the source and the green dots indicate the position of the viewers.Figure 6
Electric field distribution (Ez) for free space for a 2.45 GHz frequency. On the left, the FDTD-CPML case is shown, while on the right, the reference propagation (extended domain). In (a), (b), and (c), we show snapshots at n= 150,650,700, respectively. We observe an excellent wave absorption without presenting spurious reflected waves. The simulation lasted for 1160 time steps, and the source was turned off at n=580.
(a)
(b)
(c)The comparison of the viewer-recorded values of electric fieldEz for both the reference and the ABC cases is shown in Figure 7. To perform the analysis of the energy decay, we turned off the source at n=580 time steps (the simulation lasted for 1160 time steps). In Figure 8, we show the behaviour of such an energy within the physical domain as a function of time, demonstrating that the CPML ABC accurately dissipates the energy introduced by the source.Figure 7
Electric field behaviour in the viewers of interest, for free propagation experiment. The difference in the values of the field for the CPML and the reference solutions is negligible. When the source is turned off, the field is accurately absorbed by the CPML ABC.Figure 8
Energy decay as a function of time for free propagation experiment. It can be observed that when the source is turned off (n=580), the energy is adequately absorbed by the CPML ABC.In Table2, we present the mean squared error between the viewer values of reference and CPML solutions.Table 2
Mean squared error (Ez) for each viewer compared with the reference solution for 2.45 GHz, 5.00 GHz, and 20.00 GHz frequencies. For 5.00 GHz, we iterated up to n=3480 time steps and the source is turned off at n=1160 time steps.
Viewer
2.45 GHz
5.00 GHz
20.00 GHz
Viewer 1
3.0205 × 10–07
1.0736 × 10–07
5.0416 × 10–05
Viewer 2
5.2573 × 10–07
1.1895 × 10–07
2.2171 × 10–09
Viewer 3
4.0470 × 10–07
9.5491 × 10–08
3.0028 × 10–09


### 4.2. Parabolic Plate Reflector

In this experiment, we propagated electromagnetic waves upon the surface of a parabolic antenna, in order to obtain the distribution of electromagnetic field in such a surface as well as the reflected field in the direction of the primary radiator, which is located in the focus of the parabola.Both the physical and extended domains are of the same size than in the previous experiment. The parabolic plate is simulated to be made of silver (σ=6.3×107Ω−1m−1). The main radiator consists of a spherically symmetric source of electromagnetic field, within a hood feeder of the same material. The remaining computational domain is considered as free space. The considered geometry is shown in Figure 9, where the green dots represent the position of the numerical viewers.Figure 9
Configuration of the parabolic reflector. The parabola is described byfx=x2/0.60.4. The viewers are located in strategic points so to measure the quality of the solution. The source is located in the parabola’s focus, within the hood feeder.The electromagnetic propagation for this experiment is shown in Figure10. The behaviour of the electromagnetic field within the physical domain occurs as expected, and it features correctly the conversion from spherical to plane wavefronts. It can be clearly observed that the greater energy concentration occurs in the parabola’s focus, when it is illuminated by the plane wavefront. Figure 10 features four snapshots at different simulation times, n=100, 400, 750, and 1160, for both the CPML and the reference cases. The source was turned off at n=580 time steps.Figure 10
Distribution of electric field (Ez) for the parabolic plate experiment. At left, we observe the propagation in the CPML case, while at right, we show the reference solution. In (a), (b), (c), and (d), we present snapshots at n=100, 400, 750, and 1160 time steps. We observe an excellent absorption of waves in the CPML region, without the presence of spurious waves reflected within the physical domain. The simulation lasted n=1160 time steps, and the source was turned off at n=580.
(a)
(b)
(c)
(d)Viewers 1 and 3 are located in such positions (see Figure9) to record the field of the plane wavefront generated by the parabolic reflector. Of particular interest are viewers 5 and 6, whose position is designed to record border effects. Viewer 6 records refraction phenomena occurred in the edge of the parabolic plate, due to the energy propagated on the parabolic surface. Viewer 5 accurately records the backscattered electromagnetic field in the backwards of the parabolic plate; despite of its low intensity, this effect has already been observed [21]. Being the parabolic reflector of great interest for telecommunications engineering, its radiation pattern has been obtained in [22].Figure11 shows the behaviour of each of the viewers described above, comparing the values of the CPML and reference solutions. Figure 12 shows the energy decay after n=580 time steps (when the source is turned off), where it is possible to observe an adequate dissipation of energy within the computational domain, demonstrating the efficiency of the implemented CPML ABC.Figure 11
Behaviour of electromagnetic field in the viewers of interest, for the parabolic plate experiment. No distinguishable difference can be observed between the CPML and the reference solutions. We observe an excellent wave absorption without presenting spurious reflected waves. The simulation lasted for 1160 time steps, and the source was turned off atn=580.Figure 12
Energy decay as a function of time for the parabolic plate experiment. We observe the correct energy dissipation after the source is turned off (n=580).Table3 presents the mean square error for each viewer, between the CPML and the reference solutions.Table 3
Mean squared error (Ez) for each viewer in the parabolic reflector experiment, between the CPML and the reference solutions, for operation frequencies of 2.45 GHz, 5.00 GHz, and 20.00 GHz.
Viewer
2.45 GHz
5.00 GHz
20.00 GHz
Viewer 1
7.5934 × 10–07
3.0083 × 10–07
4.3917 × 10–08
Viewer 2
3.6161 × 10–07
1.8743 × 10–07
3.4687 × 10–08
Viewer 3
4.6658 × 10–07
2.7481 × 10–07
4.1319 × 10–08
Viewer 4
7.2940 × 10–07
5.0781 × 10–07
6.6099 × 10–08
Viewer 5
3.6715 × 10–07
3.1612 × 10–07
3.8112 × 10–08
Viewer 6
3.5750 × 10–07
1.6050 × 10–07
3.0605 × 10–08


### 4.3. Coplanar Nanowaveguide

In this section, we propagate electromagnetic waves with frequencies lying in the THz band. The basic structure of the CNWG consists of three metallic rectangular parallel plates, a central one and two at its sides, equidistantly separated from the central plate. The plates are placed on a squared dielectric substrate.Figure13 shows the configuration of the waveguide, which is characterised at a frequency of 100 THz. The three silver plates are located over a dielectric squared substrate of side 6×10−05 m, while the size of the plates is of 4.15×10−05 m long by 7.50×10−06 m width. The separation between the plates is 9.375−06 m. The orange dots represent sources and are located at a distance of 2.537×10−05 from the plates. The thickness of the substrate is 150×10−06 m. We located strategically four numerical viewers (green dots) within the physical domain, so to record electromagnetic field in such positions. Two are located between the plates to analyse the behaviour of the field within the waveguide, while the other three are located outside the waveguide. The sources were located strategically so to obtain the better propagation of the waves along the waveguide (see Figure 13). For the CNWG, the domain is extended from [6×10−05, 6×10−05] m to [6.5×10−04, 6.5×10−04] m.Figure 13
Configuration of the geometry of the coplanar nanowaveguide.Figure14 shows electric field Ez propagating through the CNWG. We observe the field propagates symmetrically along the waveguide. TE mode provides low dispersion in CNWGs, so it is optimal to build circuits and wideband electronic components. Figure 14 presents three snapshots at different time steps, for both the CPML and the reference solutions.Figure 14
Distribution of electric field (Ez) in the CNWG for a propagation frequency of 100 THz. Cases (a), (b), and (c) represent snapshots of the propagation at n=100, 350, and 1160 time steps.
(a)
(b)
(c)Figure15 shows the behaviour of field in the interest viewers, and in Table 4, we report the mean squared error for each viewer between the reference and the CPML solutions.Figure 15
Behaviour of electric field in the viewers of interest. No differences between the values of the viewers of both the reference and the CPML solutions can be observed. When turning off the source (n=580), the CPML ABC adequately absorbs arriving waves.Table 4
Mean squared error (Ez) for each viewer on the CNWG, between the reference and the CPML experiments, for a propagation frequency of 100 THz.
Viewer
100 THz
Viewer 1
3.090578701762399 × 10–08
Viewer 2
1.025905938998774 × 10–08
Viewer 3
1.025905938998774 × 10–08
Viewer 4
1.294024700943559 × 10–08
Viewer 5
1.294024700943559 × 10–08


### 4.4. Indoor Propagation

In this numerical experiment, we propagate Wi-Fi waves (ƒ=2.45GHz) in an indoor setup of two contiguous offices, based on the experiments reported in [23]. We characterised the different materials that are usually found in office indoors such as metals, wood, glass, and plastics, by their conductivity, permeability, and permittivity. Ten different materials were identified and characterised and are shown with a colour code in Figure 16, which is featured along with the corresponding values of their physical parameters in Table 5. For this experiment, the domain is of 3.728 m × 5.696 m. Due to the complexity and the size of physical domain with respect to the propagated wavelength as well as due to RAM requirements, we did not perform the reference solution.Figure 16
Media distribution within the office indoor herein studied. The values of the physical parameters characterising each material are shown in Table5.Table 5
Physical parameters of the materials considered for the indoor propagation depicted in Figure16.
Material
Conductivity
Permittivity
Permeability
σ (Ω−1m−1)
ε (Fm−1)
μ (NA−2)
Iron
9.93e6
2000ε0
0.9999μ0
Aluminium 1
37.7e6
1.0ε0
0.9991μ0
Aluminium 2
37.7e6
1.0ε0
1.00002μ0
Wood 1
10e−10
4.0ε0
0.9999μ0
Wood 2
10e−10
4.0ε0
0.999μ0
Glass
10e−12
6.0ε0
0.9991μ0
Paper 1
10e−11
7.0ε0
0.999μ0
Paper 2
10e−11
7.0ε0
0.9999μ0
Plastic 1
10e−16
2.7ε0
0.999μ0
Plastic 2
10e−16
2.7ε0
0.9999μ0Figure17 shows the electric field Ez through four snapshots at different time steps. In order to probe the quality of the CPML ABC, its stability is tested by extending the time domain up to n= 100,000 time steps and by turning off the source at n=50,000. With this numerical experiment setup, we show in Figure 18 the energy decay due to the presence of the CPML ABC, demonstrating its correct operation for large time domains.Figure 17
Distribution of electric fieldEz through the indoor setup. Panels (a), (b), (c), and (d) cases show snapshots at n=450, 1050, 2050, and 3050 time steps.
(a)
(b)
(c)
(d)Figure 18
Energy decay as a function of time for indoor propagation. We observe the accurate dissipation of energy due to the presence of the CPML ABC, after the source is turned off atn=50,000 time steps.In Figures19 and 20, we show the behaviour of electric field (Ez) recorded at selected viewers in both offices, respectively. As it can be clearly observed, after the source is turned off, the field clearly decreases as a function of time, fact that confirms the stability of the ABC herein implemented, for both long time domains and a wide range of propagating materials. As it can be observed in such figures, the six viewers show a very interesting feature; notwithstanding its positions on the indoor setup, for early times in the simulation, they show the transients before arriving to a quasistationary field propagation status. Such transients are also present when the source is turned off at n= 50,000 time steps, arriving to another quasistationary state of no field, after the energy is dissipated by the CPML ABC.Figure 19
Behaviour of the electric fieldEz in the viewers of the first office.Figure 20
Behaviour of the electric fieldEz in the viewers of the second office.Viewers in the lower office (viewers 1, 2, and 3) exhibit greater values in the recorded electric field values, mainly due to the fact that the source is located within this office, so the field is more intense within this part of the propagation domain. The first viewer shows that it receives the first wavefront coming directly from the source, which travels and leaves the viewer number one. But after certain time, the viewer records the arriving of further wavefronts coming from some edges that scattered waves towards the viewer 1 direction, as well as some reflected waves entering into this region.As viewer 2 is located in a free space region, adjacent to two boundary regions with a CPML ABC, the amount of scattered and reflected waves arriving at this viewer is very much lower, which is clearly shown in the fact that the field recorded by this viewer is more regular than the field recorded by viewers 1 and 3. This is not the case for viewer 3, which is closer to several edges and material interfaces, so it is subject to record wavefronts at different time steps, coming from different scatter and reflection processes.For the case of viewers 4, 5, and 6, located in the upper (second) office, the behaviour of the recorded field is different. These viewers clearly record lower values of the electric fieldEz, as the waves arrive attenuated to them. Particularly, viewer 5 behaves as viewer 2, provided that it is surrounded by two CPML ABCs, which prevents most of scattered and reflected waves to arrive to it. Viewers 4 and 6 show clearly different wavefronts arriving from different processes of reflection and scattering of waves from the surroundings.


## 4.1. Free Space Propagation

The physical domain set up is 1.25 m × 2.25 m while the reference (extended) domains are of 20 m × 20 m, 10 m × 10 m, and 5 m × 5 m for 20.00 GHz, 5.00 GHz, and 2.45 GHz, respectively, with the source located in the centre. Numerical viewer locations are shown in Figure5. The results of the simulation for the CPML and reference solutions are shown in Figure 6, where the absorption in the CPML region is observed.Figure 5
Configuration for the free space propagation scenario. Orange lines indicate the beginning of the absorption region (CPML); the orange dot in the centre indicates the location of the source and the green dots indicate the position of the viewers.Figure 6
Electric field distribution (Ez) for free space for a 2.45 GHz frequency. On the left, the FDTD-CPML case is shown, while on the right, the reference propagation (extended domain). In (a), (b), and (c), we show snapshots at n= 150,650,700, respectively. We observe an excellent wave absorption without presenting spurious reflected waves. The simulation lasted for 1160 time steps, and the source was turned off at n=580.
(a)
(b)
(c)The comparison of the viewer-recorded values of electric fieldEz for both the reference and the ABC cases is shown in Figure 7. To perform the analysis of the energy decay, we turned off the source at n=580 time steps (the simulation lasted for 1160 time steps). In Figure 8, we show the behaviour of such an energy within the physical domain as a function of time, demonstrating that the CPML ABC accurately dissipates the energy introduced by the source.Figure 7
Electric field behaviour in the viewers of interest, for free propagation experiment. The difference in the values of the field for the CPML and the reference solutions is negligible. When the source is turned off, the field is accurately absorbed by the CPML ABC.Figure 8
Energy decay as a function of time for free propagation experiment. It can be observed that when the source is turned off (n=580), the energy is adequately absorbed by the CPML ABC.In Table2, we present the mean squared error between the viewer values of reference and CPML solutions.Table 2
Mean squared error (Ez) for each viewer compared with the reference solution for 2.45 GHz, 5.00 GHz, and 20.00 GHz frequencies. For 5.00 GHz, we iterated up to n=3480 time steps and the source is turned off at n=1160 time steps.
Viewer
2.45 GHz
5.00 GHz
20.00 GHz
Viewer 1
3.0205 × 10–07
1.0736 × 10–07
5.0416 × 10–05
Viewer 2
5.2573 × 10–07
1.1895 × 10–07
2.2171 × 10–09
Viewer 3
4.0470 × 10–07
9.5491 × 10–08
3.0028 × 10–09


## 4.2. Parabolic Plate Reflector

In this experiment, we propagated electromagnetic waves upon the surface of a parabolic antenna, in order to obtain the distribution of electromagnetic field in such a surface as well as the reflected field in the direction of the primary radiator, which is located in the focus of the parabola.Both the physical and extended domains are of the same size than in the previous experiment. The parabolic plate is simulated to be made of silver (σ=6.3×107Ω−1m−1). The main radiator consists of a spherically symmetric source of electromagnetic field, within a hood feeder of the same material. The remaining computational domain is considered as free space. The considered geometry is shown in Figure 9, where the green dots represent the position of the numerical viewers.Figure 9
Configuration of the parabolic reflector. The parabola is described byfx=x2/0.60.4. The viewers are located in strategic points so to measure the quality of the solution. The source is located in the parabola’s focus, within the hood feeder.The electromagnetic propagation for this experiment is shown in Figure10. The behaviour of the electromagnetic field within the physical domain occurs as expected, and it features correctly the conversion from spherical to plane wavefronts. It can be clearly observed that the greater energy concentration occurs in the parabola’s focus, when it is illuminated by the plane wavefront. Figure 10 features four snapshots at different simulation times, n=100, 400, 750, and 1160, for both the CPML and the reference cases. The source was turned off at n=580 time steps.Figure 10
Distribution of electric field (Ez) for the parabolic plate experiment. At left, we observe the propagation in the CPML case, while at right, we show the reference solution. In (a), (b), (c), and (d), we present snapshots at n=100, 400, 750, and 1160 time steps. We observe an excellent absorption of waves in the CPML region, without the presence of spurious waves reflected within the physical domain. The simulation lasted n=1160 time steps, and the source was turned off at n=580.
(a)
(b)
(c)
(d)Viewers 1 and 3 are located in such positions (see Figure9) to record the field of the plane wavefront generated by the parabolic reflector. Of particular interest are viewers 5 and 6, whose position is designed to record border effects. Viewer 6 records refraction phenomena occurred in the edge of the parabolic plate, due to the energy propagated on the parabolic surface. Viewer 5 accurately records the backscattered electromagnetic field in the backwards of the parabolic plate; despite of its low intensity, this effect has already been observed [21]. Being the parabolic reflector of great interest for telecommunications engineering, its radiation pattern has been obtained in [22].Figure11 shows the behaviour of each of the viewers described above, comparing the values of the CPML and reference solutions. Figure 12 shows the energy decay after n=580 time steps (when the source is turned off), where it is possible to observe an adequate dissipation of energy within the computational domain, demonstrating the efficiency of the implemented CPML ABC.Figure 11
Behaviour of electromagnetic field in the viewers of interest, for the parabolic plate experiment. No distinguishable difference can be observed between the CPML and the reference solutions. We observe an excellent wave absorption without presenting spurious reflected waves. The simulation lasted for 1160 time steps, and the source was turned off atn=580.Figure 12
Energy decay as a function of time for the parabolic plate experiment. We observe the correct energy dissipation after the source is turned off (n=580).Table3 presents the mean square error for each viewer, between the CPML and the reference solutions.Table 3
Mean squared error (Ez) for each viewer in the parabolic reflector experiment, between the CPML and the reference solutions, for operation frequencies of 2.45 GHz, 5.00 GHz, and 20.00 GHz.
Viewer
2.45 GHz
5.00 GHz
20.00 GHz
Viewer 1
7.5934 × 10–07
3.0083 × 10–07
4.3917 × 10–08
Viewer 2
3.6161 × 10–07
1.8743 × 10–07
3.4687 × 10–08
Viewer 3
4.6658 × 10–07
2.7481 × 10–07
4.1319 × 10–08
Viewer 4
7.2940 × 10–07
5.0781 × 10–07
6.6099 × 10–08
Viewer 5
3.6715 × 10–07
3.1612 × 10–07
3.8112 × 10–08
Viewer 6
3.5750 × 10–07
1.6050 × 10–07
3.0605 × 10–08


## 4.3. Coplanar Nanowaveguide

In this section, we propagate electromagnetic waves with frequencies lying in the THz band. The basic structure of the CNWG consists of three metallic rectangular parallel plates, a central one and two at its sides, equidistantly separated from the central plate. The plates are placed on a squared dielectric substrate.Figure13 shows the configuration of the waveguide, which is characterised at a frequency of 100 THz. The three silver plates are located over a dielectric squared substrate of side 6×10−05 m, while the size of the plates is of 4.15×10−05 m long by 7.50×10−06 m width. The separation between the plates is 9.375−06 m. The orange dots represent sources and are located at a distance of 2.537×10−05 from the plates. The thickness of the substrate is 150×10−06 m. We located strategically four numerical viewers (green dots) within the physical domain, so to record electromagnetic field in such positions. Two are located between the plates to analyse the behaviour of the field within the waveguide, while the other three are located outside the waveguide. The sources were located strategically so to obtain the better propagation of the waves along the waveguide (see Figure 13). For the CNWG, the domain is extended from [6×10−05, 6×10−05] m to [6.5×10−04, 6.5×10−04] m.Figure 13
Configuration of the geometry of the coplanar nanowaveguide.Figure14 shows electric field Ez propagating through the CNWG. We observe the field propagates symmetrically along the waveguide. TE mode provides low dispersion in CNWGs, so it is optimal to build circuits and wideband electronic components. Figure 14 presents three snapshots at different time steps, for both the CPML and the reference solutions.Figure 14
Distribution of electric field (Ez) in the CNWG for a propagation frequency of 100 THz. Cases (a), (b), and (c) represent snapshots of the propagation at n=100, 350, and 1160 time steps.
(a)
(b)
(c)Figure15 shows the behaviour of field in the interest viewers, and in Table 4, we report the mean squared error for each viewer between the reference and the CPML solutions.Figure 15
Behaviour of electric field in the viewers of interest. No differences between the values of the viewers of both the reference and the CPML solutions can be observed. When turning off the source (n=580), the CPML ABC adequately absorbs arriving waves.Table 4
Mean squared error (Ez) for each viewer on the CNWG, between the reference and the CPML experiments, for a propagation frequency of 100 THz.
Viewer
100 THz
Viewer 1
3.090578701762399 × 10–08
Viewer 2
1.025905938998774 × 10–08
Viewer 3
1.025905938998774 × 10–08
Viewer 4
1.294024700943559 × 10–08
Viewer 5
1.294024700943559 × 10–08


## 4.4. Indoor Propagation

In this numerical experiment, we propagate Wi-Fi waves (ƒ=2.45GHz) in an indoor setup of two contiguous offices, based on the experiments reported in [23]. We characterised the different materials that are usually found in office indoors such as metals, wood, glass, and plastics, by their conductivity, permeability, and permittivity. Ten different materials were identified and characterised and are shown with a colour code in Figure 16, which is featured along with the corresponding values of their physical parameters in Table 5. For this experiment, the domain is of 3.728 m × 5.696 m. Due to the complexity and the size of physical domain with respect to the propagated wavelength as well as due to RAM requirements, we did not perform the reference solution.Figure 16
Media distribution within the office indoor herein studied. The values of the physical parameters characterising each material are shown in Table5.Table 5
Physical parameters of the materials considered for the indoor propagation depicted in Figure16.
Material
Conductivity
Permittivity
Permeability
σ (Ω−1m−1)
ε (Fm−1)
μ (NA−2)
Iron
9.93e6
2000ε0
0.9999μ0
Aluminium 1
37.7e6
1.0ε0
0.9991μ0
Aluminium 2
37.7e6
1.0ε0
1.00002μ0
Wood 1
10e−10
4.0ε0
0.9999μ0
Wood 2
10e−10
4.0ε0
0.999μ0
Glass
10e−12
6.0ε0
0.9991μ0
Paper 1
10e−11
7.0ε0
0.999μ0
Paper 2
10e−11
7.0ε0
0.9999μ0
Plastic 1
10e−16
2.7ε0
0.999μ0
Plastic 2
10e−16
2.7ε0
0.9999μ0Figure17 shows the electric field Ez through four snapshots at different time steps. In order to probe the quality of the CPML ABC, its stability is tested by extending the time domain up to n= 100,000 time steps and by turning off the source at n=50,000. With this numerical experiment setup, we show in Figure 18 the energy decay due to the presence of the CPML ABC, demonstrating its correct operation for large time domains.Figure 17
Distribution of electric fieldEz through the indoor setup. Panels (a), (b), (c), and (d) cases show snapshots at n=450, 1050, 2050, and 3050 time steps.
(a)
(b)
(c)
(d)Figure 18
Energy decay as a function of time for indoor propagation. We observe the accurate dissipation of energy due to the presence of the CPML ABC, after the source is turned off atn=50,000 time steps.In Figures19 and 20, we show the behaviour of electric field (Ez) recorded at selected viewers in both offices, respectively. As it can be clearly observed, after the source is turned off, the field clearly decreases as a function of time, fact that confirms the stability of the ABC herein implemented, for both long time domains and a wide range of propagating materials. As it can be observed in such figures, the six viewers show a very interesting feature; notwithstanding its positions on the indoor setup, for early times in the simulation, they show the transients before arriving to a quasistationary field propagation status. Such transients are also present when the source is turned off at n= 50,000 time steps, arriving to another quasistationary state of no field, after the energy is dissipated by the CPML ABC.Figure 19
Behaviour of the electric fieldEz in the viewers of the first office.Figure 20
Behaviour of the electric fieldEz in the viewers of the second office.Viewers in the lower office (viewers 1, 2, and 3) exhibit greater values in the recorded electric field values, mainly due to the fact that the source is located within this office, so the field is more intense within this part of the propagation domain. The first viewer shows that it receives the first wavefront coming directly from the source, which travels and leaves the viewer number one. But after certain time, the viewer records the arriving of further wavefronts coming from some edges that scattered waves towards the viewer 1 direction, as well as some reflected waves entering into this region.As viewer 2 is located in a free space region, adjacent to two boundary regions with a CPML ABC, the amount of scattered and reflected waves arriving at this viewer is very much lower, which is clearly shown in the fact that the field recorded by this viewer is more regular than the field recorded by viewers 1 and 3. This is not the case for viewer 3, which is closer to several edges and material interfaces, so it is subject to record wavefronts at different time steps, coming from different scatter and reflection processes.For the case of viewers 4, 5, and 6, located in the upper (second) office, the behaviour of the recorded field is different. These viewers clearly record lower values of the electric fieldEz, as the waves arrive attenuated to them. Particularly, viewer 5 behaves as viewer 2, provided that it is surrounded by two CPML ABCs, which prevents most of scattered and reflected waves to arrive to it. Viewers 4 and 6 show clearly different wavefronts arriving from different processes of reflection and scattering of waves from the surroundings.


## 5. Performance Tests

To test the performance of the CPML ABC herein implemented, we programmed a conventional version (algorithm shown in Figure3) and our memory-optimised version (algorithm shown in Figure 4), running the experiment with single- and double-precision variables, as well as comparing different compilers.Usage of single-precision variables yields to saving almost 50% of RAM memory as well as about 30% of processing time. Nevertheless, it might introduce numerical differences of about10−7 [24].We compared the memory saving of our memory-optimised version against the conventional one. For instance, for the parabolic reflector experiment for a propagation frequency of 2.45 GHz, the memory-optimised version requires 14.4 MB and 7.3 MB for double and single precision, respectively. The not-optimised version uses 16.7 MB and 8.7 MB in double and single precision, respectively. The version of the algorithm that uses splitted ABC (Berenger’s PML [1]) used 17.6 MB and 9.0 MB in double and single precision, respectively. We can observe a memory saving of 14% with respect to the nonoptimised version, as well as an 18% with respect to the nonsplit version, for double precision. For single precision, we save 50% for all the cases. Nevertheless, the more the discrete points, the memory saving is about 30% with respect to nonsplit version and about 40% with respect to the split version.Table6 shows the RAM memory allocated by the more memory-demanding experiments. As it can be observed, there is a memory saving of 18.4%, 18.3%, and 15.6% with double-precision variables for free space, parabolic reflector, and CNWG, respectively, while for the single-precision mode, the memory savings are of about 18.4%, 18.3%, and 15.8%, respectively.Table 6
Amount of allocated memory used by the optimised version in double and single precision, for the different numerical experiments. The used memory is reported by the operating system (Linux kernel 3.13).
Numerical experiment
Double-C
Single-C
Double-E
Single-E
Empty space (20 GHz)
+Optimised
691.6 MB
353.2 MB
3.7 GB
1.9 GB
∗Nonoptimised
846.7 MB
431.1 MB
4.5 GB
2.3 GB
Parabolic plate (20 GHz)
+Optimised
691.8 MB
353.6 MB
15.9 GB
8.0 GB
∗Nonoptimised
844.8 MB
419.9 MB
19.7 GB
9.8 GB
Waveguide (100 THz)
+Optimised
27.3 MB
14.2 MB
2.4 GB
1.2 GB
∗Nonoptimised
32.3 MB
16.7 MB
2.9 GB
1.5 GBTherefore, the saving memory version saves at least 15% of RAM memory. This percentage increases as the computational domain increases too. Single precision offers more memory saving, although its precision limits to seven digits, so it could introduce numerical errors. In Figure21, we show the absolute error produced by the usage of single-precision variables with respect to the double-precision ones; the parabolic reflector, the CNWG, and the indoor propagation numerical experiments absolute errors are featured in Figures 21(a), 21(b), and 21(c), respectively.Figure 21
Error induced by single-precision variables with respect to double-precision ones.
(a)
Parabolic reflector numerical experiment.L2 norm error is of 7.403291 × 10−07 (2.45 GHz, snapshot at n=1160 time steps (7.09 s))
(b)
CNWG numerical experiment.L2 norm error is of 2.630922 × 10−06 (100 THz, snapshot at n=1161 time steps (1.740295218690000 × 10−04 s))
(c)
Indoor propagation numerical experiment.L2 norm error is of 3.0817 × 10−11 (100 THz, snapshot at n=100,000 time steps). It can be clearly observed that for such long time periods, single-precision version is numerically stableTo test computing time in serial mode, we performed experiments on a quad-core Intel Core i5-4670 at 3.40 GHz (without HT), over a Linux kernel 3.13.0-65-generic SMP operating system. Comparing memory-optimised versus conventional versions, there is no observable difference in computing times. Nevertheless, some differences between different compilers are observed. For the parabolic reflector at 20 GHz, by compiling with gfortran 4.9.2 with optimisation and vectorisation flags activated, we obtain computing times of 10 m, 55 s and 04 m, 36 s for the memory-optimised version in double and single precision, respectively. By compiling with PGI FORTRAN, we obtain computing times of 11 m, 02 s and 05 m, 13 s for the conventional version for double and single precision. The difference in computing times between the compilers is negligible; actually, memory-optimised version generates more cycles, providing a similar performance.Table7 features computing times using the memory-optimised version for double and single precision, for the different numerical experiments herein performed, where we can observe a computing time saving of about 40%.Table 7
Computing time required in double- and single-precision modes, for the different numerical experiments, considering CPML and reference solutions.
Numerical experiment
Double-C
Single-C
Double-E
Single-E
Free space (20 GHz)
+Optimised
04 m, 27 s
01 m, 58 s
23 m, 11 s
10 m, 31 s
∗Nonoptimised
04 m, 41 s
02 m, 16 s
25 m, 14 s
12 m, 11 s
Parabolic plate (20 GHz)
+Optimised
10 m, 27 s
04 m, 34 s
04 h, 0 m, 29 s
01 h, 45 m, 05 s
∗Nonoptimised
10 m, 54 s
04 m, 36 s
04 h, 0 m, 29 s
01 h, 41 m, 29 s
Waveguide (100 THz)
+Optimised
28 s
10 s
12 m, 24 s
04 m, 26 s
∗Nonoptimised
29 s
11 s
14 m, 35 s
05 m, 32 s


## 6. Conclusions

We analysed the behaviour and performance of CPML ABC version by Komatitsch and Martin [3], for different frequencies in the band from MHz to THz, observing an excellent absorption and energy dissipation and lacking spurious reflected waves being propagated within the physical domain. Moreover, we found ad hoc CPML parameters that work perfectly for all numerical experiments herein performed, for the whole range of explored frequencies. It ought to also be remarked that this CPML ABC implementation features stability for long time periods. Our implementation of the FDTD-CPML algorithm provides at least a 15% of memory saving versus a conventional implementation, without sacrificing performance. We also find that single-precision mode saves about 50% of RAM memory and computing time, without introducing sensitive numerical error due to the usage of single-precision variables.Furthermore, although for the sake of clarity, our implementation of the herein developed algorithm was done in 2D, it can be straightforwardly extended to 3D just by building the vectorskx, ax, and bx in (13), (14), and (15) carefully in each boundary of the problem for them to correspond to the absorbed physical field on each boundary, taking into account their corresponding physical parameters. It ought to be remarked that as our implementation of the CPML ABC saves memory through the drastic reduction of the allocated arrays with respect to conventional PML ABC implementations, the achieved optimisation in 3D could be of at least 64% [18].



---

*Source: 5710943-2018-04-02.xml*