International Journal of Antennas and Propagation
Volume 2018 (2018), Article ID 5710943, 20 pages
https://doi.org/10.1155/2018/5710943

Research Article

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

A. Rodríguez-Sánchez,¹ C. Couder-Castañeda,² J. J. Hernández-Gómez,² I. Medina,² S. Peña-Ruiz,³ J. Sosa-Pedroza,³ and M. A. Enciso-Aguilar³

¹Universidad del Valle de Mexico, Campus Hispano, vía José López Portillo No. 346 y 352, San Lorenzo Tetlixtac, 55700 Coacalco de Berriozabal, Mexico
²Instituto Politécnico Nacional, Centro de Desarrollo Aeroespacial, Belisario Domínguez 22, Col. Centro, Del. Cuauhtémoc, 06010 Ciudad de México, Mexico
³Instituto Politécnico Nacional, Escuela Superior de Ingeniería Mecánica y Eléctrica, Av. Luis Enrique Erro S/N, Unidad Profesional Adolfo López Mateos, Zacatenco, Delegación Gustavo A. Madero, 07738 Ciudad de México, Mexico

Correspondence should be addressed to C. Couder-Castañeda; ccouder@ipn.mx

Received 14 August 2017; Revised 8 December 2017; Accepted 25 December 2017; Published 2 April 2018

Academic Editor: Renato Cicchetti

Copyright © 2018 A. Rodríguez-Sánchez 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.

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.

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, most solvers 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 Section 2, 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]

where and 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 . The external current and electric charge densities, and , respectively, are the sources of electromagnetic field, and they can be functions of space and time, that is, and .

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 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

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 , then the finite differences version of such equations are, respectively,

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 , with , 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 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 ( or ) field in the direction for each time step as

For instance, the auxiliary equation for takes the form

The vectors , , and are estimated, within the absorption region, as

where = , is the CPML thickness, , and is the number of layers of the CPML ABC. It is important to note that within the physical domain, , , and ; within the CPML

where is the vacuum permittivity,

is the damping function profile,

is the shift in the frequency domain [5], and and are the degree of polynomial variation for and , respectively.

The vectors , , and are calculated in the intermediate steps , while the vectors , , and are calculated in the intermediate steps .

Equations (8), (9), and (10) with the implemented CPML take the form

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 propagation

To 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 Table 1. 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 and ((14) and (15)), which in turn fully depends on , and ((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 , , and , which from (13), (16), and (17), turn to be , , , , and . 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 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 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

3.0

1.0

3.0

0.08

20

The 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 the direction is IMAX and in the direction is JMAX, then the whole domain including the CPML might be labelled in the direction with and in the direction with , where , , , , and 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 Figure 3, 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 as , , , and .

Figure 4: Flux diagram of the numerical algorithm using memory variables only in the absorption region. The cycles are numbered as .

The spatial increments and are estimated considering ( being the propagated wavelength), while is estimated with the Courant stability criterion . The sources of electromagnetic field are modelled as , where is the frequency and 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 energy within the physical domain, at time step , is given by

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 Figure 5. 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.

(a)

(b)

(c)

Figure 6: Electric field distribution () 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 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 .

The comparison of the viewer-recorded values of electric field 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 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 (), the energy is adequately absorbed by the CPML ABC.

In Table 2, we present the mean squared error between the viewer values of reference and CPML solutions.

Table 2: Mean squared error () 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 time steps and the source is turned off at 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 (). 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 by . 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 Figure 10. 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, 100, 400, 750, and 1160, for both the CPML and the reference cases. The source was turned off at time steps.

(a)

(b)

(c)

(d)

Figure 10: Distribution of electric field () 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 , 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 time steps, and the source was turned off at .

Viewers 1 and 3 are located in such positions (see Figure 9) 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].

Figure 11 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 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 at .

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 ().

Table 3 presents the mean square error for each viewer, between the CPML and the reference solutions.

Table 3: Mean squared error () 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.

Figure 13 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 m, while the size of the plates is of m long by m width. The separation between the plates is m. The orange dots represent sources and are located at a distance of from the plates. The thickness of the substrate is 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 [, ] m to [, ] m.

Figure 13: Configuration of the geometry of the coplanar nanowaveguide.

Figure 14 shows electric field 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.

(a)

(b)

(c)

Figure 14: Distribution of electric field () in the CNWG for a propagation frequency of 100 THz. Cases (a), (b), and (c) represent snapshots of the propagation at , 350, and 1160 time steps.

Figure 15 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 (), the CPML ABC adequately absorbs arriving waves.

Table 4: Mean squared error () 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 () 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 Table 5.

Table 5: Physical parameters of the materials considered for the indoor propagation depicted in Figure 16.

Material
Conductivity
Permittivity
Permeability

()
()
()

Iron
9.93e6
2000
0.9999

Aluminium 1
37.7e6
1.0
0.9991

Aluminium 2
37.7e6
1.0
1.00002

Wood 1
10e−10
4.0
0.9999

Wood 2
10e−10
4.0
0.999

Glass
10e−12
6.0
0.9991

Paper 1
10e−11
7.0
0.999

Paper 2
10e−11
7.0
0.9999

Plastic 1
10e−16
2.7
0.999

Plastic 2
10e−16
2.7
0.9999

Figure 17 shows the electric field 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 100,000 time steps and by turning off the source at . 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.

(a)

(b)

(c)

(d)

Figure 17: Distribution of electric field through the indoor setup. Panels (a), (b), (c), and (d) cases show snapshots at , 1050, 2050, and 3050 time steps.

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 at time steps.

In Figures 19 and 20, we show the behaviour of electric field () 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 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 field in the viewers of the first office.

Figure 20: Behaviour of the electric field 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 field , 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 Figure 3) 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 about [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.

Table 6 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 GB

Therefore, 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 Figure 21, 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.

(a) Parabolic reflector numerical experiment. L² norm error is of 7.403291 × 10⁻⁰⁷ (2.45 GHz, snapshot at time steps (7.09 s))

(b) CNWG numerical experiment. L² norm error is of 2.630922 × 10⁻⁰⁶ (100 THz, snapshot at time steps (1.740295218690000 × 10⁻⁰⁴ s))

(c) Indoor propagation numerical experiment. L² norm error is of 3.0817 × 10⁻¹¹ (100 THz, snapshot at time steps). It can be clearly observed that for such long time periods, single-precision version is numerically stable

Figure 21: Error induced by single-precision variables with respect to double-precision ones.

To 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.

Table 7 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 vectors , , and 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].

Conflicts of Interest

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

Acknowledgments

The authors acknowledge EDI grant by SIP/IPN. This research was partially supported by projects SIP 20160105, 20171536, 20170721, and 20171027 and IPN. The authors also acknowledge FS0005 computational time grant at FSF by CDA/IPN.

References

J.-P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” Journal of Computational Physics, vol. 114, no. 2, pp. 185–200, 1994.

J. A. Roden and S. D. Gedney, “Convolution PML (CPML): an efficient FDTD implementation of the CFS–PML for arbitrary media,” Microwave and Optical Technology Letters, vol. 27, no. 5, pp. 334–339, 2000.

D. Komatitsch and R. Martin, “An unsplit convolutional perfectly matched layer improved at grazing incidence for the seismic wave equation,” Geophysics, vol. 72, no. 5, pp. SM155–SM167, 2007.

R. Martin, D. Komatitsch, and A. Ezziani, “An unsplit convolutional perfectly matched layer improved at grazing incidence for seismic wave propagation in poroelastic media,” Geophysics, vol. 73, no. 4, pp. T51–T61, 2008.

R. Martin and C. Couder-Castañeda, “An improved unsplit and convolutional perfectly matched layer absorbing technique for the navier-stokes equations using cut-off frequency shift,” Computer Modeling in Engineering and Sciences, vol. 63, no. 1, pp. 47–77, 2010.

R. Luebbers, F. P. Hunsberger, K. S. Kunz, R. B. Standler, and M. Schneider, “A frequency-dependent finite-difference time-domain formulation for dispersive materials,” IEEE Transactions on Electromagnetic Compatibility, vol. 32, no. 3, pp. 222–227, 1990.

P. M. Goorjian and A. Taflove, “Direct time integration of Maxwell’s equations in nonlinear dispersive media for propagation and scattering of femtosecond electromagnetic solitons,” Optics Letters, vol. 17, no. 3, pp. 180–182, 1992.

R. Vuduc and K. Czechowski, “What GPU computing means for high-end systems,” System, vol. 8, p. 10, 2011.

T. Zhang, Y. Du, T. Huang, and X. Li, “GPU-accelerated 3D reconstruction of porous media using multiple-point statistics,” Computational Geosciences, vol. 19, no. 1, pp. 79–98, 2014.

J. Owens, D. Luebke, N. Govindaraju et al., “A survey of general-purpose computation on graphics hardware,” Computer Graphics Forum, vol. 26, no. 1, pp. 80–113, 2007.

N. Nakata, T. Tsuji, and T. Matsuoka, “Acceleration of computation speed for elastic wave simulation using a graphic processing unit,” Exploration Geophysics, vol. 42, no. 1, pp. 98–104, 2011.

J. Gomes, G. Teodoro, A. Melo, J. Kong, T. Kurc, and J. Saltz, “Efficient irregular wavefront propagation algorithms on Intel(R) Xeon Phi(TM),” in 2015 27th International Symposium on Computer Architecture and High Performance Computing (SBAC-PAD), pp. 25–32, Florianopolis, Brazil, 2016.

G. Civario, S. Delaney, and M. Lysaght, “Preparing a seismic imaging code for the Intel Knights Landing Xeon Phi processor,” Advances in Parallel Computing, vol. 27, pp. 585–590, 2016.

A. Heinecke, A. Breuer, M. Bader, and P. Dubey, “High order seismic simulations on the Intel Xeon Phi processor (Knights Landing),” Lecture Notes in Computer Science, vol. 9697, pp. 343–362, 2016.

M. Hernández, J. Cebrián, J. Cecilia, and J. García, “Evaluation of 3-d stencil codes on the Intel Xeon Phi coprocessor,” Advances in Parallel Computing, vol. 27, pp. 197–206, 2016.

F. Mori, M. Matsumoto, and T. Furumura, “Performance optimization of the 3D FDM simulation of seismic wave propagation on the Intel Xeon Phi coprocessor using the ppOpen-APPL/FDM library,” Lecture Notes in Computer Science, vol. 8969, pp. 66–76, 2015.

J. D. Jackson, Classical Electrodynamics, John Wiley & Sons. Inc., New York, NY, USA, 1999.

R. Martin and D. Komatitsch, “An unsplit convolutional perfectly matched layer technique improved at grazing incidence for the viscoelastic wave equation,” Geophysical Journal International, vol. 179, no. 1, pp. 333–344, 2009.

M. Arroyo, C. Couder-Castañeda, A. Trujillo-Alcantara, I.-E. Herrera-Diaz, and N. Vera-Chavez, “A performance study of a dual Xeon-Phi cluster for the forward modelling of gravitational fields,” Scientific Programming, vol. 2015, Article ID 316012, 14 pages, 2015.

C. Couder-Castañeda, H. Barrios-Piña, I. Gitler, and M. Arroyo, “Performance of a code migration for the simulation of supersonic ejector flow to SMP, MIC, and GPU using OpenMP, OpenMP+LEO, and OpenACC directives,” Scientific Programming, vol. 2015, Article ID 739107, 20 pages, 2015.

E. Jull, “Aperture antennas and diffraction theory,” in Electromagnetics and Radar Series, Peter Peregrinus Ltd, Stevenage, UK, 1981.

J. Sosa Pedroza, M. Enciso Aguilar, S. Peña Ruiz, A. Rodríguez Sánchez, and E. Garduño Nolasco, “A new proposal to convertnear-field into far-field for parabolic reflectors,” Microwave and Optical Technology Letters, vol. 58, no. 7, pp. 1595–1599, 2016.

J. Sosa, S. Coss, A. Rodríguez et al., “Indoor 2.4 GHz microwave propagation study using 3D FDTD approach,” Electronics Letters, vol. 47, no. 24, pp. 1308-1309, 2011.

C. Couder-Castañeda, C. Ortiz-Alemán, M. Orozco-Del-Castillo, and M. Nava-Flores, “TESLA GPUs versus MPI with OpenMP for the forward modeling of gravity and gravity gradient of large prisms ensemble,” Journal of Applied Mathematics, vol. 2013, Article ID 437357, 15 pages, 2013.