Block algorithm for the joint difference solution of the d'Alembert and Maxwell's equations

Abstract

A characteristic feature of mathematical modeling at the present stage of development is the consideration of the architecture of the computer system, not only for stage of compiling a computer program, but also during the development of a numerical method and synthesis of the mathematical model. This method significantly broadens the researcher's ability to search for the optimal mapping of the numerical method to the mentioned architecture, in the sense of accelerating computations. In this paper, this idea is illustrating by examples of the basic mathematical model of computational electrodynamics and optics, Maxwell's equations, and the FDTD. This modification allows to reducing the data exchange rate between the operational and cache memory due to the greater number of arithmetic operations per one grid function in solving the d'Alembert equation. On the other hand, freely use the technologies FDTD method and ready-made software implementations for setting the incident wave, imposing the absorbing layers, taking into account the dispersion of the medium.

1. Introduction

Despite the old beginning of researches [ 1 ] in the field of the numerical solution of the equations of Maxwell, interest in this subject domain only increases over time. This remark first of all belongs to application of a method of final differences [ 2 ] – to version of the numerical decision most popular so far the specified equations. In the seventies the last century [ 3 ] he has received own abbreviation of FDTD (Finite-Difference Time-Domain) under which it is widely known to this day.

The relevance of development of a FDTD method is due to several reasons. First, broad demand in new subject domains: to nanophotonics [ 4 ], radiobiology [ 5 ], etc. owing to community of mathematical model. Indeed, by means of Maxwell's equations are described any processes connected with the wave nature of electromagnetic radiation without restrictions. Secondly, the necessity of taking into account the features of modern computing architectures: cluster [ 6 ], vector [ 7 ], etc. If previously the performance of computers was increased mainly due to the increase in the clock speed of the central processor, but now it is increasing using various methods of parallel processing of data. In the study proposed by the authors of this work, the emphasis in synthesizing a new variant of the FDTD method is made on the basis of the hierarchical structure of computer storage devices, in particular, on the possibility of optimizing communications between the operational and cache memory of the processor. Specified subject sparingly illuminated in the scientific press in force the interdisciplinary nature of the task. Experts in calculus mathematics traditionally are interested in other problems of the theory of differential schemes: increase in an order of approximation of a differential task [ 10 ], designing of mobile net areas [ 11 ]. In turn, developers of mathematical software that implements the FDTD method, seeking to keep pace with the development of "large" modern forms of parallel and distributed computing: using processors [ 12 ] the clouds [ 13 ].

Nevertheless, fresh works appear on the block algorithms for the difference solution of the basic equalizations of electrodynamics [ 14-16 ]. This testifies to the attention to "small" forms of organization of calculations, taking into account the inconspicuous architectural features of the processors and demonstrates significant acceleration at the expense of such an account. In computational practice, it has long been possible, with the example of matrix calculations, to manage the duration of calculations by changing the block sizes [ 17, 18 ] in block algorithms. Unfortunately, the features of the theory of difference schemes (the absence of the need for multiplication of dense matrices, the traditional emphasis on stability problems) until recently hampering the penetration of blockiness into this subject area.

Another obstacle to the development of block methods is the complexity of controlling the loading/unloading of data into the processor's memory cache during the computing process. Load (but not unload) statements are present only in Assembly language and their use is not prescriptive: data may not be loaded into the cache after they are executed if the block is already there.In more common programming languages, the management of communications between the operational and cache memory is achieved indirectly by receiving a partitioning of the cyclic constructs referred to as tiling [ 19, 20 ].

The application of this technique for the joint difference solution of the equations of d’Alembert and Maxwell and is devoted to the proposed publication.

2. Peculiarities of the joint difference solution of the d'Alembert and Maxwell’s equations

The theory of the joint difference solution of the d'Alembert and Maxwell’s equations is described in detail in [ 21, 22 ], the results of his experimental research, significant for the chosen subject. The following software and hardware tools were used during the experiments: the Intel Core i7-3770 processor, the Ubuntu 16.04.1 operating system (the 4.4 kernel), the gcc 5.3 compilerand Meep 1.3 (compiled for comparison by the same compiler) is a free and open-source software package for simulating electromagnetic systems via the finite-difference time-domain (FDTD) method. Actually being a reference for a wide range of researchers [ 23 ]. The grid area was chosen with the dimensions of 10000 × 10000 nodes in space and 200 nodes in time, which provided sufficient memory load (as for real computational experiments in nanophotonics) for an acceptable duration of calculations.

In the case of simulating the propagation of TM-mode (in terms of the fundamental work [ 24 ]) differential solution of the equations of d'alembert and Maxwell differ from each other on the magnitude of the machine precision and converge to the analytical. When using the MEEP package, the duration of the calculations made 124.75 sec., calculations on the author's software implementation of the difference solution of Maxwell's equations lasted 112.71 sec., for the difference solution of the d’Alembert equation lasted 41.34 sec. Acceleration in the case of the difference solution of the d’Alembertequation in 41.34 sec. and 2.73 times, respectively, can not be explained by a decrease in computational complexity by 10% in the difference solution of the d’Alembert equations compared to that for Maxwell. The authors associate the observed effect with an increase in the ratio of the number of arithmetic operations that fall on the calculation of the differential pattern, to the amount of memory involved in working with the same pattern. For a two-dimensional Yee scheme, the specified value is 3/2 (three operations are performed on two different field projections); for the first and second difference equations and 7/4 (7 operations on 3 field projections and the value of dielectric permeability) for the third; for the second scheme – 9/3 (9 operations on two different time layers of the same projection and the value of dielectric permeability), which apparently leads to a significant decrease in the intensity of communications between the operational and cache memory.

The joint difference solution combines the advantages of both approaches: acceleration for the dalamber equation and the developed Toolkit (overlays of absorbing layers, the setting of the incident wave, etc.) for the Maxwell’s equations. The Maxwell’s equations were experimentally solved in PML layers [ 24 ] with a thickness of 100 knots along the edges of the grid area, and the dalamber equation was solved in the center. The duration of the transactions in this case amounted to 47.47 sec. and slightly superior to the previous result (in which the calculations in PML layers was not carried out), obtained by accelerations of 2.63 times compared to the Meep package and 2.37 times compared to the author's implementation of the Yee scheme.

Further, the authors use the observed effect of reducing the duration of calculations with a decrease in the intensity of communication between the operational and cache memory for further acceleration of calculations in the case of a joint solution.

3. Block algorithm for the joint difference solution

To date, block algorithms for the difference solution of Maxwell's equations [ 14 ] (2009) and for d’Alembert [ 16 ] (2015) are known. The method of transition to cloud computing with diamond toroidal block shape (Diamond Torre Algorithm) proposed in the last work served as the basis for the synthesis of the author's block algorithm of joint difference solution. The transition to cloud computing consists of the following two steps.

The first one performs calculations on theh time layers of the grid area to update the values of the grid functions in the sub-area of the PML location. In figure 1 it can be said "bucket handle". Due to the information dependence on the iteration space coinciding with the grid area, the values in the nodes adjacent to the distance h and smaller to the PML, which are already related to the wave equation, are also to be calculated at this stage.So in the node separated from the absorbing layers on k other nodes ( k ≤ h ) the value of the grid function on the ktime layer will be formed.All such units will make up the "bucket wall" in figure 1.

The second stage is characterized by the organization of calculations in the internal volume of the "bucket" according to the algorithm, which the authors called the wave by analogy with the method of transition to the block of [ 19 ]. As shown in figure 2, during the computing process inside the bucket, the values of the grid functions are calculated from left to right in such a way that the front of the process has the form of an inclined plane.In the nodes of the grid area before it, the function values are not yet defined, in the nodes of the grid area after it are found on the layer h. On the very surface of the front, values are calculated on different layers in ascending order with increasing height and decreasing node abscissa. After the second stage of the algorithm, the first one comes again, and their alternation continues until the values of the grid functions on all layers of the region are found.

In contrast to the diamond toroidal algorithm from [ 16 ], such an approach is easier to implement and leads to a slight but stable reduction in the duration of calculations compared to the other method of transition to blockage. T1 T2

V

It is noteworthy that the best results for both algorithms are achieved for h=20the maximum wave height (or torus) at which the block still fits into the cache memory as a whole (the volume of L3 cache for Intel Core i7-3770 is 8 MB). Indeed, ath>20 the block is no longer placed in the cache and you have to load it in parts, and for h<20 the fast cache memory is not used entirely. Both of these circumstances lead to an increase in communications and, as a consequence, the total duration of calculations.

4. The study of the silicon Fresnel’s lens

Of particular interest is the calculation and creation of elements to control the radiation of the Novosibirsk free electron laser (NovoFEL) [ 25 ], which is the most powerful source of terahertz radiation. By now, focusing elements in the form of binary diffraction lenses made of silicon [ 26, 27 ] are known for it, which are characterized by the known advantages and disadvantages of binary optics. Diffraction elements favorably with refractive [ 28 ] much smaller thickness, high radiation resistance, greater manufacturability in use and manufacture. However, the focusing element with the binary relief mentioned in [ 26, 27 ] is characterized by low efficiency (often not exceeding 40% [ 28 ]), in which most of the energy of the illuminating beam does not fall into the focus area.

However, diffraction elements with a continuous profile, free from this drawback, have been known for a long time. Moreover, the calculation of the binary lens itself is traditionally [ 28 ] preceded by the calculation of the Fresnel’s lens, which is such an element. The problem of forming a continuous relief on a silicon plate has not yet been solved technologically, the processes of mechanical and chemical processing of silicon are considered to be sufficiently coarse to obtain an uninterrupted profile with the necessary accuracy. However, advances in plasma etching of the diamond surface (not inferior to silicon in hardness and inertness to various kinds of effects) [ 29 ] allow us to hope for the rapid improvement of the technology of reactive-ion etching of silicon demonstrated so far in the manufacture of binary relief [ 30 ].

The model, numerical method and software complex presented in [ 22, 31, 32 ] and based on the joint difference solution of Maxwell's equations with the use of block algorithm of organization of calculations are chosen as a tool for calculating diffraction on various silicon lenses. Вased on [ 26,27,30 ], we assume the refractive index of silicon n=3.42 for the wavelength λ=141 µm, to which we further give all distances. Thus, the aperture of the refractive lens put equal 101.74 λ, the thickness and the radius of curvature of 10.33λ and 130.41λ, respectively. Then, according to the geometric optics [ 33 ], the focus area will be located at a distance f=50.87λ from the right pole. Calculation of Fresnel’s lens traditionally assumed the height of its relief equal to h=λ/(n-1), which is unacceptable in this case. For a lens with such a high numerical aperture, half of its Fresnel’s zones will be less than the width of the wavelength, which is why the work of the element is strongly unpredictable. Taking h=2.07λ, Fresnel's own lens was calculated (figure3), consisting of 5 zones, the extreme of which is characterized by an acceptable width of 5.18 λ.

Figure 3 shows the lens focusing of the radiation incident from left to right. On figure 4, there is a noticeable drop in the intensity I (normalized to the intensity of the incident beam) for the lens under discussion as compared to the refractive one. The intensity after the first one is δ=0.9 from the intensity value after the second one in the coordinates of the maxima. The latter are shifted by 4.73 λ and δ=2.47 λ to the left for refractive and diffraction lenses, respectively.

Of particular interest is the modeling of elements with manufacturing process errors [ 34 ], especially since such studies have not yet been conducted for silicon elements. Let's assume that the height error is associated with overheating of the plate surface (figure5), leading to an uneven increase in the height of the profile. On figure 5, the value of the discussed value is taken to be equal to λ for better clarity, then consider the more real smaller values (table 2). Indeed, the maximum height error for the spent process is unlikely to exceed 10%.

The results of modeling the Fresnel’s lenses with manufacturing process errors presented in table 2 were somewhat unexpected. The natural assumption about the negative influence of errors on the profile height on the focusing efficiency was justified only for sufficiently large error values exceeding a quarter of the wavelength. The deviation of 0.15λ (figure 4, table 2) led to an increase in intensity in the focus area by a quarter compared to the refractive lens. This may be due to the high Fresnel reflection coefficient of silicon, and, consequently, to the significant fluctuations of the intensity of the transmitted wave with a slight change in the height of the profile. Curiously, for a plane-parallel plate, the maximum intensity also fell by a width of 0.15λ.

References

[1] Cron G 1944 Equivalent circuit of the field equations of Maxwell Proc . IRE 32 289 - 299

[2] Yee K S 1966 Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media IEEE Trans . Antennas Propag 14 302 - 307

[3] Taflove A and Brodwin M 1975 Numerical solution of steady-state electromagnetic scattering problems using the time-dependent Maxwell's equations IEEE Transactions of Microwave Theory and Techniques 2 ( 8 ) 623 - 630

[4] Gavrilov A V , Golovashkin D L , Doskolovich L L , Dyachenko P N , Khonina S N , Kotlyar V V , Kovalev A A , Nalimov A G , Nesterenko D V , Pavelyev V S , Shuyupova Y O , Skidanov R V and Soifer V A 2014 Diffracrive Nanophotonics (Boca Raton: CRC Press) p 679

[5] Perov S Yu and Bogacheva E V 2014 Theoretical and experimental dosimetry in the evaluation of the biological action of electromagnetic fields of wearable radio stations . Flat phantoms Radiation biology: Radioecology 54 ( 1 ) 57 - 61

[6] Guiffaut C and Mahdjoubi K 2001 A parallel FDTD algorithm using the MPI library IEEE Antennas and Propagation Magazine 43 ( 2 ) 94 - 103

[7] Vorotnikova D G and Golovashkin D L 2014 CUBLAS-aided Long Vector Algorithms Journal of Mathematical Modelling and Algorithms in Operations Research 13 ( 4 ) 425 - 431

[8] Malysheva S A and Golovashkin D L 2016 Implementation of the FDTD algorithm on GPU using a pyramid method Computer Optics 40 ( 2 ) 179 - 187 DOI: 10.18287/ 2412 -6179-2016- 40-2- 179 -187

[9] Vorotnikova D G and Golovashkin D L 2017 Difference solutions of the wave equation on GPU with reuse of pairwise sums of the differential template Computer Optics 41 ( 1 ) 134 - 138 DOI: 10.18287/ 2412 -6179-2017-41-1- 134 -138

[10] Xiao F 2007 High-order accurate split-step FDTD method for solution of Maxwell's equations Electronics Letters 43 ( 2 ) 72

[11] Fidel B , Heyman E , Kastner R and Zioklowski R W 1997 Hybrid ray-FDTD moving window approach to pulsepropagation Journal of Computational Physics 138 ( 2 ) 480 - 500

[12] FDTD solver (Access mode : http://www.acceleware.com/fdtd-solvers)

[13] EMA3D Version 3 . 3 parallel modification (Access mode: http://www .ema3d.com)

[14] Orozco D and Guang G 2009 Mapping the FDTD Application to Many-Core Chip Architectures Parallel Processing 309-316

[15] Grosser T and Cohen A 2013 Split Tiling for GPUs: Automatic Parallelization Using Trapezoidal Tiles Proc . of the GPGPU 6 24 - 31

[16] Perepelkina A Yu and Levchenko V D 2015 Diamond Torre Algorithm for High-Performance Wave Modeling Keldysh Institute Preprints 18 20

[17] Golub G H 1996 Matrix Computations (Baltimor: Johns Hopkins University Press) p 726

[18] Demmel J W 1997 Applied Numerical Linear Algebra (Baltimor : Johns Hopkins University Press) p 435

[19] Wolfe M 1986 Loops skewing: The wavefront method revisited International Journal of Parallel Programming 15 ( 4 ) 279 - 293

[20] Wolfe M 1989 More Iteration Space Tiling Proceedings of the ACM 655-664

[21] Golovashkin D L and Yablokova L V 2012 Joint finite-difference solution of the d'Alembert and Maxwell's equations . One-dimensional case Computer Optics 36 ( 4 ) 527 - 533

[22] Golovashkin D L , Buldygin E Y and Yablokova L V 2014 Joint finite - difference solution of the d'Alembert and Maxwell's equations . Two-dimensional case Computer Optics 381 20 - 27

[23] Oskooi A F , Roundyb D and Ibanescua M 2010 Meep: A flexible free-software package for electromagnetic simulations by the FDTD method Computer Physics Communications 181 687 - 702

[24] Taflove A and Hagness S 2005 Computational Electrodynamics: The Finite-Difference TimeDomain Method (Boston: Arthech House Publishers) p 1006

[25] Knyazev B A , Kulipanov G N and Vinokurov N A 2010 Novosibirsk terahertz free electron laser: instrumentation development and experimental achievements Measurement Science and Technology 21 13

[26] Agafonov A N , Volodkin B O , Kaveev A K , Knyazev B A , Kropotov G I , Pavelyev V S , Soyfer V A , Tukmakov K N , Tsygankova E V and Choporova Yu Yu 2013 Optoelectronics, Instrumentation and Data Processing 49 ( 2 ) 98 - 105

[27] Agafonov A N , Vlasenko M G , Volodkin B O , Gerasimov V V , Kaveev A K , Knyazev B A , Kropotov G I , Pavelyev V S , Palchikova I G , Soyfer V A , Stupak M F , Tukmakov K N , Tsygankova E V and Choporova Yu Yu 2013 Diffraction lenses for powerful bundles of terahertz radiation News of RAS Series physical 77(9 ) 1360 - 1362

[28] Golovashkin D L , Doskolovich L L , Kazansky N L , Kotlyar V V , Pavelyev V S , Skidanov R V , Soyfer V A and Khonina S N 2007 Diffraction computer optics (Moscow: Fizmatlit) p 736

[29] Karlsson M and Nikolajeff F 2003 Diamond micro -optics: microlenses and antireflection structured surfaces for the infrared spectral region Optics Express 11 ( 5 ) 502 - 507

[30] Agafonov A N , Volodkin B O , Kaveev A K , Knyazev B A , Kropotov G I , Pavelyev V S , Tukmakov K N and Choporova Yu Yu 2014 Control of a terahertz laser radiation of cross and modal composition by means of elements of binary silicon optics Computer Optics 38 ( 4 ) 763 - 769

[31] Yablokova L V and Golovashkin D L 2018 Block algorithms for the joint difference solution of the D'Alembert and Maxwell equations Computer Optics 42 ( 2 ) 320 - 327 DOI: 10.18287/ 2412 - 6179-2018-42-2- 320 -327

[32] Certificate of state registration of the computer program № 2017613903 Joint difference solution of the equations of D'Alembert and Maxwell

[33] Butikov E I 2003Optics (St. Petersburg : Nevsky Dialect) p 480

[34] Golovashkikn D L , Volkov A V , Eropolov V A , Kazanskiy N L , Karpeev S V , Moiseev O Y , Pavelyev V S , Artushenko V G and Kashin V V 2007 Studying fabrication errors or the diffraction grating on the end face of a silver-halide fiber Optical Memory & Neural Networks Information Optics 16 ( 4 ) 263 - 268