Application of the pyramid method in difference solution d'Alembert equations on graphic processor with the use of Matlab

Abstract

The paper proposes a modification of the pyramid method for constructing algorithms for the difference solution of the d'Alembert equation on a graphics processor in the event of a shortage of video memory. The authors demonstrate the effectiveness of the method on the practical example of dividing the grid area into two sub domains. Acceleration reaches the characteristic for the case of a domain entirely located in the video memory. In the article investigated the effectiveness of using the author's approach depending on the height of the pyramid and showed the boundaries of applicability of the proposed modification.

Keywords

• The method of the pyramids
• the grid area
• difference solution
• computing acceleration

1. Introduction

E k 1  2Eik  E k 1 i i ht2  c2 Eik1  2hEz2ik  Eik1 ( 1 ) is written with respect to the grid function defined on the domain Dh  {(tk , zi ) : tk  kht , k  0,1,..., Nt  T ht , zi  iht , i  1,..., N z  Lz strength is, c is the speed of light in free space, T and Lz are size of the region in time and space. h  1} , where E the value of the electric field z

High-Performance Computing / L.V. Yablokova, D.L. Golovashkin Below is a fragment of the computational procedure for solving ( 1 ) in MATLAB, where c1  c2 h2 h2 , c5  2 c ht  . t z % Placement of grid functions on two time layers in video memory

E1=zeros(1,Nz,'gpuArray'); E2=zeros(1,Nz,'gpuArray'); for k=1:2:Nt % Passage through time layers of the grid area through one E1(2:Nz-1)=2*E2(2:Nz-1)-E1(2:Nz-1)+c1*diff(E2,2); % Calculations on the layer k

E1( 2 )=sin(c5*k); % Hard radiation condition on the layer к E2(2:Nz-1)=2*E1(2:Nz-1)-E2(2:Nz-1)+c1*diff(E1,2); % Calculations on the layer k+1

E2( 2 )=sin(c5*(k+1)); % Hard radiation condition on the layer k+1

End E=gather(E2); % Transfer of results to RAM

For Nz  5107 and Nt  100 the duration of calculations on the Intel Core i7 CPU was 57.08 s., On the GeForce GTX TITAN X GPU - 5.55 s. (acceleration of 10.29 times) using MATLAB 2015b and the operating system CentOS 7.2. Both used arrays occupied 762 MB in memory, however, during the computations on the CPU, the memory requirements increased by one and a half time, on the GPU the memory requirements increased threefold. Apparently, with the implementation of calculations for the design E1(2:Nz-1)=2*E2(2:Nz-1)-E1(2:Nz-1)+c1*diff(E2,2) on the CPU, the execution of the operation of numerical differentiation diff(E2,2) resulted in allocating additional memory for two copies of the array E2, and the execution on the GPU of the design as a completely required separate area of memory for all the arrays involved and for double copying E2. MATLAB takes about 0.4 gigabytes in RAM, but does not use video memory for its placement. Thus, the execution of the whole algorithm on the CPU was accompanying by the extraction of 1.52 GB. In addition, the execution of the whole algorithm on the GPU was accompanying by the extraction of 2.24 GB. Moreover, if the researcher has a video card with 2 GB of memory (like most popular video processors now) then the organization of calculations on the GPU becomes impossible. In his previous publication [7], using the difference scheme for the Maxwell equations, the CUDA C software tool the authors proposed to solve this problem using the method of pyramids.

3. The pyramid method application

2 .

Will this be possible in this case, given that MATLAB is not specialized for working with graphics processors and its tools in this area are very meager?

The essence of the mentioned method in the author's modification consists in splitting the grid domain into overlapping sub regions that fit in the video memory completely, with the subsequent organization of communications in the production of vector computations in each sub region separately. In this case, transfers from RAM to video and vice versa are performed not on each time layer, but through a certain number of them h (the height of the pyramid). This, on the one hand, leads to a reduction in h the number of communications. On the other hand, to the duplication of arithmetic operations in overlapping fragments of grid subdomains (the form of pyramids is available in the two-dimensional case).

A fragment of the computational procedure implementing this strategy in the case under consideration is presented below, where N  Nz % creating temporary layers on CPU and GPU

E1=zeros(1,Nz); E2=zeros(1,Nz); E1m=zeros(1,h); E2m=zeros(1,h); E1g=zeros(1,N+h,'gpuArray'); E2=zeros(1,N+h,'gpuArray');

for t=1:h:Nt % Passage through the pyramids % work with the left subdomain

E1g=gpuArray(E1(1:N+h)); E2g=gpuArray(E2(1:N+h)); % Forwarding CPU ==> GPU for k=1:2:h % Calculations inside the first pyramids

E1g(2:N+h-k)=2*E2g(2:N+h-k)-E1g(2:N+h-k)+c1*diff(E2g(1:N+h-k+1),2);

E1g( 2 )=sin(c5*(t+k-1));

High-Performance Computing / L.V. Yablokova, D.L. Golovashkin E2g(2:N+h-k-1)=2*E1g(2:N+h-k-1)-E2g(2:N+h-k-1)+c1*diff(E1g(1:N+h-k),2);

E2g( 2 )=sin(c5*(t+k)); end E1(2:N-h)=gather(E1g(2:N-h)); E2(2:N-h)=gather(E2g(2:N-h)); % Forwarding GPU ==> CPU

E1m(1:h)=gather(E1g(N-h+1:N)); E2m(1:h)=gather(E2g(N-h+1:N)); % work with the right subdomain

E1g(1:N+h-1)=gpuArray(E1(N-h+1:Nz)); E2g(1:N+h-1)=gpuArray(E2(N-h+1:Nz)); for t=1:2:h % Calculation by layers of the pyramid

E1g(t+1:N+h-2)=2*E2g(t+1:N+h-2)-E1g(t+1:N+h-2)+c1*diff(E2g(t:N+h-1),2);

E2g(t+2:N+h-2)=2*E1g(t+2:N+h-2)-E2g(t+2:N+h-2)+c1*diff(E1g(t+1:N+h-1),2); end E1(N+1:Nz-1)=gather(E1g(h+1:N+h-2)); % Forwarding GPU ==> CPU E2(N+1:Nz-1)=gather(E2g(h+1:N+h-2));

E1(N-h+1:N)=E1m(1:h); E2(N-h+1:N)=E2m(1:h); % Replenishment of the result end

In the course of experiments with the new algorithm, the dependence of the calculation time on the height of the pyramid. The Table 1 contains the results.

Height of the pyramid, h

Computation time (s)

Acceleration 2 4 10 20 50 53.02 29.58 15.49 10.75 7.9 1.08 1.93 3.7 5.31 7.23

4. Conclusion

Acknowledgements References

Thus, the method of pyramids can be effectively using in arranging calculations for solving difference equations with the help of MATLAB on graphic processors in the case when arrays storing values of grid functions do not fit into video memory as a whole. The development of the proposed algorithm for cases of large dimensions will be the next stage of the authors' research in this direction.

The research leading to these results has received funding from the Russian Science Foundation grant №16-47-630560-r_a.

References

[1] Taflove A , Hagness S. Computational Electrodynamics: The Finite-Difference Time-Domain Method . Boston: Aerotech House Publishers, 2005 ; 1006 p.

[2] Elsherbeni А , Demir V . The Finite-Difference Time-Domain Method for Electromagnetics with MATLAB Simulations . Scitech Publishing Inc. 2009 ; 426 p.

[3] Grigorjev IS , Mejlihov EZ . Physical Values: Reference Book. Moscow: Energoatomizdat Publisher, 1991 ; 1232 p. (in Russian) .

[4] Malysheva SA , Golovashkin DL . Realization of the difference solution of the Maxwell equations on graphic processors by the pyramid method . Computer Optics 2016 ; 40 ( 2 ): 179 - 187 . DOI: 10 .18287/ 2412 -6179-2016-40-2- 179 -187.

[5] Kozlova ES , Kotlyar VV . Simulation of the propagation of a short two-dimensional pulse of light . Computer Optics 2012 ; 36 ( 2 ): 158 - 164 .

[6] Vorotnikova DG , Golovashkin DL . Difference solution of the wave equation on graphical processors with repeated use of pairwise sums of the differential pattern . Computer Optics 2017 ; 41 ( 1 ): 134 - 138 . DOI: 10 .18287/ 2412 -6179-2017-41-1- 134 -138.

[7] Vorotnikova DG , Golovashkin DL . Algorithms with "long" vectors for solving grid equations of explicit difference schemes . Computer Optics 2015 ; 39 ( 1 ): 87 - 93 . DOI: 10 .18287/ 0134 -2452-2015-39-1- 87 -93.