-This paper proposes an innovative algorithm for running in a distributed and parallel fashion the well-known Finite-Difference Time-Domain method. Given the dependence among data, for the proposed distributed version, care has been taken in the approach splitting the involved data. This is to avoid a great deal of data exchange among different hosts, so as to minimize the communication overhead. Performance comparisons between the proposed solution and the sequential one has shown that improvements can be achieved using the proposed version running on several hosts. Moreover, the distributed version let overcome the problem of central memory exhaustion, typical of this method, due to the large matrices to be handled. Index Terms-FDTD, parallelization, distribution, algorithms, data dependence
The FDTD Finite-Difference Time-Domain method [
The original formulation consists of a sequential algorithm
executing a possibly huge number of iterations on a single core
of a host. Distribution is highly desirable for solving the typical
problems on the domain under investigations, i.e.
electromagnetic fields, etc., due to the high amount of memory used and
computation necessary. However, a distributed version of the
algorithm has to minimize communication overhead due to
data exchange among hosts. When the distributed version is
available, then an underlying platform can be put to use for
achieving real parallelism and distribution [
A proper analysis of the data dependencies among the
several parts of the initial sequential algorithms is therefore
needed to achieve a distributed version minimising data
exchanges. In this field, several automatic or semi-automatic
techniques can assist in finding structural and data
dependencies [
Copyright c 2017 held by the authors.
This paper analyses the most important aspect of distributing the typical FDTD algorithm and shows the performance improvements when the distributed version is used. The following sections of the paper describe, firstly the mathematical background, secondly the sequential version of the algorithm, thirdly an initial optimization for the algorithm, fourthly the distributed and parallel version. Finally, conclusions are drawn.
The technique called FIT (Finite Integration Theory [
These equations can be integrated on dual meshes at interlaced time intervals, based on the leap-frog scheme as shown in Fig. 1.
This method leads to the following equations, representing an iterative scheme.
D~ xn+ 21 (k) B~yn+1(k + ) = ~ E
D~ xn 21 (k) 1 ~
D "r B~yn(k + 1 ~
B r t c0 [Hyn(k + 1 z 2 )
Hyn(k )] t c0 [E~xn+ 21 (k + 1) z
E~xn+ 12 (k)] Fig. 1. x
Discretization in one-dimension with the dual grids scheme. S = z, S~ = y z, where x and y are the discretization steps.
The above iterative equations are updated (time marching) as shown in the leap-frog scheme of Fig. 2. The steps marked in Fig. 2 with 2 and 4 are relative to the costitutive relations of the fields and are punctual operations. While the steps marked with 1 and 3 are relative to the Maxwell’s curl equations (the data used for these updates are located in adjacent cells).
Of course, an electromagnetic wave can not propagate at a speed higher than the light velocity. Therefore, the time step must not exceed the time that the wave would take to propagate between 2 points on the grid. In 2D the wave propagation direction changing the traveling distance must be taken into account (see Fig. 3 and 4).
In 3D in order to implement the iterative scheme we use the Yee cell that is shown in the Fig. 4.
A. The Gaussian Pulse
fg(t) = e ( twdgg )2
Where the dg component represents the time delay, while the wg parameter generically is the width of the pulse. The corresponding peak value is obtained for t = dg. If we express the dg delay and the wg pulse width, depending on the time step t, we have the discrete version of the previous equation: (1)
In the Gaussian pulse equation in the discrete version (2), the t temporal time is not included explicitly.
The time-discrete version of the equation is appropriate to be translated into an algorithm which can be computed for a given set of data. Finally, it is necessary to take into account other features of the physical problem, that is the following three constraints.
1) The electric boundary condition, since in the internal points of a perfect conductor the electric field E vanishes. 2) The magnetic boundary condition, which can ba easily handled if the computational domain stops at magnetical nodes. 3) The radiation condition, for the ability to absorb waves without reflection. pulse = exp(-0.5*(pow((t0-T)/2.0))); dz[IC][JC] = pulse;
The study of the physical phenomenon and its mathematical functions lead us to outline its implementation. The algorithm takes as input several fundamental values for the study of this phenomenon, including 5 matrices, initialized to the needed values for the computation, and initialized in an optimized way (see the following). Once values are received, the algorithm executes some checks and then reaches the 4 fundamental for loops, and the pulse variable calculation.
The inputs for the algorithm are: 1) the initial matrices ga[][], dz[][], ez[][], hx[][] and hy[][], containing the necessary values required to perform the calculation; 2) the number of steps nsteps to be performed, expressing the number of times that the algorithm has to be executed; 3) T the absolute time, from which the study of the physical phenomenon starts; 4) the source parameter t0, essential for calculating the pulse variable (usually set to 20.0); 5) the spread of the source, needed for calculating the pulse variable (set to 6.0).
There are two control loops, an external while and an inner for. The first external loop checks whether or not the condition that ends the execution of the algorithm (value 0 or a negative number) has been reached. It is a condition taken as an input.
The second loop checks the number of nsteps that have to be executed by the algorithm. It also responsible for increasing the T absolute time variable, which measures the time from the zero instant of the beginning of the experiment until to the instant n where the experiment ends, and also is used to compute the electromagnetic values. The time variable is a float, since time flows with steps of 0.5.
The subsequent 4 for loops compute 4 of the 5 starting matrices (ga[][] is constant). The first loop computes dz[][] matrix, using hy[][] and hx[][] matrices, multiplied by the constant 0.5. With respect to the elements of matrix hy[][] to be computed, the previous row of matrix hy[][] and the next column of hx[][] have to be made available. The dependence from such data has great implications for the transformation of the algorithm into its distributed version.
Then, element by element there will be some additions to compute the dz[][] matrix, as in the following: dz[i][j] = dz[i][j] + 0.5*(hy[i][j] hy[i-1][j] - hx[i][j] + hx[i][j-1]);
Before the second for loop, a block of instructions sets the pulse variable according to time T. After that, the resulting value is assigned to a given portion of dz[][], previously calculated. This is said the Gaussian impulse of the source. The corresponding code is as follows:
The second for loop calculates matrix ez[][], by means of the matrices ga[][] and dz[][], (the second matrix is updated by the first for loop and the gaussian impulse). The code is as follows: ez[i][j] = ga[i][j] * dz[i][j];
The third for loop calculates the matrix hx[][], adding it to the multiplication between the constant value 0.5 and the updated ez[][] matrix. The next column of the matrix is required (with further implications, as in the previous cases). The relative code is: hx[i][j] = hx[i][j] + 0.5 * (ez[i][j] ez[i][j+1]);
Finally, the fourth for loop calculates the matrix hy[][], like in the third loop, by using matrix hy[][], the constant value 0.5 and the updated matrix ez[][]. The corresponding code is: hy[i][j] = hy[i][j] + 0.5 * (ez[i+1][j] - ez[i][j]);
The said algorithm has been developed in C++, and even the sequential version has better performances than some corresponding commercial software systems. The simulation times for our developed tool were significantly lower than those required by simulations–using the same hardware. This reduction is extremely advantageous when complex electromagnetic structures are studied.
Firstly, it is required to initialize the 5 matrices that represent the electromagnetic magnitudes under study. They are: ga[][], dz[][], ez[][], hx[][], hy[][]. Both the size of the matrices and the values they have to have, will be given as input to the algorithm. In particular, the size of all matrices must be identical between them, while the values can be different. This derives from the mathematical function used to study the physical phenomenon. Such values can be given in a configuration file, in order to make the program more versatile and suitable for the investigation of the phenomenon from different starting conditions.
The dynamic initialization of the 5 starting matrices is provided through the initMatrix function, which has as input the matrix to initialize, its size, and the initial value that the matrix elements must take. These elements are float types. initMatrix function consists of 2 for loops, one for loop scans the rows and one for loop scans the columns. Within the 2 for loops, it is performed an indexed assignment of the matrices to the desired value. It can be expected that a small part of the array will be set to a starting value other than the rest (which is usually 0 or 1.0).
A performance comparison between the proposed dynamic initialization and a static initialization of the matrices is shown in table II. Compilation and execution average times are given for both versions.
If the matrix size increases, a slight improvement in performances of the static initialized matrices has been observed. This result is apparently positive, considering the big size typical of the matrices, however it only affects about 0.7% of the initialization time, and would depend on the matrix size. In addition, this approach has several limits, some of which are really considerable. A brief list of such limits is reported below: 1) The compilation time, which must be actually added to the runtime program for each run, is in percentage really high compared to the run time (due to the huge amount of data needed for pre-storage). This consideration, combined with the fact that data are not constant, but taken as input, will be such that this time has to be considered whenever the program is executed. 2) The limit on the number of elements that can be compiled, which causes a virtual memory consumption error (depending on the host, and for a standard host being around 853 elements). 3) The dependence on a c++ compiler (or g++) that each host must have installed on itself to let the program work. 4) The inability to create an executable version of the software system, which in fact reduces the user-friendliness for less experienced users.
To further improve the dynamic initialization performance, we have used pointers to access matrix elements, since the arithmetic pointer is faster than array indexing. k = rows * columns; for (i = 0; i <k; i++) *(matrix + i) = value;
The matrix variable (a float pointer) holds the onedimensional array representing the matrix; the integers i and j are used to scan the matrix; the row and column are int variables holding the total number of rows and columns; and finally the value contains the float value of the cells in the array, given to initMatrix function. Further tests indicated an improvement in performance, hence reducing the slight gap between the two approaches.
By distributing the 5 starting matrices on multiple hosts it is
possible to parallelize and distribute execution, overcome the
physical limits of a single host (large amount of data in volatile
memory, etc.), and speed up the processing [
The intrinsic structure of the implementation (see the above description of the sequential version), limits the distribution approach. In order to compute matrices (i.e. dz[][], hx[][] and hy[][]) some rows or columns are needed from, previous or successive, other matrices (hy[][], hx[][] and ez[][]). The parallelism of the problem at hand can be seen according to the producer-consumer model. Once the row or column values have been calculated, which are needed by another host to calculate other sub-matrices, these are sent by means of sockets. The underlying idea for distribution is to have either a chain of n hosts, or a number of hosts organised as a grid. Moreover, each host running a part of the processing manages the communication with the nearest hosts.
The first part of the execution is handled by server.cpp program that chooses the number of hosts that will be used (among the available n), according to the size of the 5 starting matrices: ga[][], dz[][], ez[][], hx[][] and hy[][]. Once the x hosts have been chosen, the above matrices are split into x sub-matrices, and sent to the x hosts. The sub-matrices are spread depending on the position of the host.
The criteria for splitting the 5 matrices (ga[][], dz[][], ez[][], hx[][] and hy[][]) into x submatrices are mainly two: the first concerns the x quantity of the submatrices to be created, the second the arrangement of them.
The quantity x of the sub-matrices, which is equal to the x number of the n computers to be available, depends on two factors: the size of the starting matrices, and a function that determines the accurate “multiple” of computers and their logical interconnection.
According to some size ranges of the matrices, memory occupation and execution times on a sample host are observer, and then a number y is advised as a possible way to handle the data. The size split will provide an exact multiple of hosts (see Fig. 6, in which each box represents both a submatrix and a host).
This number of hosts and parts is described by the following two mathematical equations
y + 1 y + 1 y2 + 2y + 1 x1 = =
2 2 4 where y is the number suggested according to the initial occupancy test on a sample host. Moreover, the following equation will be used.
y y + 2 y2 + 2y x2 = =
2 2 4
The equation pair gets the y value and returns the number of the sub-matrices in which the starting matrices have to be split. The number given by the first equation x1 represents the number of rows where submatrices will be arranged, whereas the second equation, giving x2 defines the number of columns. To make the system more efficient and increase parallelism, it is necessary to reduce drastically the waiting times of all computers. By arranging submatrices in this way, each host that finishes its calculations, will provide inputs to, and activate, both the host to its right (right), and the underlying host (lower), according to the scheme shown in Fig. 7.
This system significantly reduces the times for the global calculation of starting matrices compared to the chain method, as shown in the following example reported in Fig. 8.
At the time t1 (for instance equal to 10 minutes), the first host finishes its calculations, and passes the results to the right
Fig. 8. Execution times using the chain system and using the grid approach. Although the first, calculating the same amount of data and using the same number of hosts it takes 40 minutes total, the second 30. and the bottom host, ending their processing at the instant t2 (20 minutes). Both of them finally relocate their results to the last host, which ends up to calculate their data at t3 time (equal to 30 minutes).
Increasing the number of sub-matrices in which the starting matrices are divided, the times for the grid version of the algorithm improve, in contrast with the time for the chain version.
This paper described the development of a distributed system for calculating 5 electromagnetic magnitudes ((ga[][], dz[][], ez[][], hx[][] and hy[][]). It was necessary to initialize the matrices to well-defined values, through the dynamic initialization. It was then through the server.cpp program that the 5 starting matrices have been split into as many sub-matrices as the selected x hosts among the available n. Then data are transferred to the client.cpp program, and then are synchronized. The sockets offered extreme flexibility and ease of use for handling communication. At this point, the necessary computation has been spread among several hosts, each having a components performing the needed task.