Numerical Simulation of Heat Transfer in Semiconductor Heterostructures Karine K. Abgaryan1,2, Ilya S. Kolbin1 1 Dorodnitsyn Computing Center FRC CSC RAS, Moscow, Russia 2 Moscow Aviation Institute (NRU) (MAI), Moscow, Russia, kristal83@mail.ru Abstract The paper considers the problems of building-up numerical models of heat transfer in multilayer nanostructure AlAs/GaAs. The calculations relied that the temperature inside each layer is equal through its full width, while the changes proceed in discrete steps at the interfaces. The problem solving was carried out on the basis of hybrid finite-difference meshless method, for time discretization an implicit finite-difference scheme was chosen. In these settings dimensional approximation was performed by applying the radial basis function package. The choice of weight coefficients for meshless solution was carried on through the system of linear equations, and the centers and the “widths” of the bases were fixed. The obtained results were compared for different bases (Gaussian function, multiquadric and inverse multiquadric) together with and without normalization approximating model. It was shown that normalization gives higher accuracy results in some cases. For the quality assessment of the obtained results, they were compared with those, calculated according classical method of finite differences. The offered algorithm gave similar results. Investigation the matter of increasing speed method using concurrent computation revealed that application of multiprocessor system can allow to reach significant growth of efficiency of the developed algorithm. 1 Introduction Nowadays there has been a significant amount of interest in multilayer heterostructures with nanoscale layer thickness. In a point of fact, these structures have become a frequent practice in micro and nanoscale semiconductor electronics production. The reason for this is high values of thermal conduction, electrical conductivity etc. as compared to traditional materials. Thus, simulation of procedures of physical processes in corresponding structures is an urgent science and technology task. One of the promising approaches for heat transfer simulation is the usage of meshless methods. Application of corresponding algorithms allows to penetrate some difficulties commonly found in classical methods of finite differences, finite elements, and etc. such as exponential increase of complexity considering the problems of high dimensionality, solving problems with compound range and etc. [3]. The widely accepted meshless approach is approximation of simulating processes by means of radial basic functions (RBF) [19]. Several equivalent forms of recording approximating RBF. This paper presents the following: u ( x ) =  i =1  i  i ( x ) , M (1)  i ( x ) =  ( ri ( x ),  i ) , (2) ri ( x ) = x − x ic , (3) whereω – weight coefficients, ε – “widths” of bases, xc – centers coordinates. For example, Hardy multiquadrics are used as radial functions φ [7, 8, 17]: Copyright © 2019 for the individual papers by the papers' authors. Copyright © 2019 for the volume as a collection by its editors. This volume and its papers are published under the Creative Commons License Attribution 4.0 International (CC BY 4.0). In: Sergey I. Smagin, Alexander A. Zatsarinnyy (eds.): V International Conference Information Technologies and High- Performance Computing (ITHPC-2019), Khabarovsk, Russia, 16-19 Sep, 2019, published at http://ceur-ws.org 37 Computer Design of New Materials ___________________________________________________________________________________________  ( x) = 1 +  2 r 2 . (4) The algorithm, suggested by Kanza [9, 11, 14, 16], involves selection of weight coefficients of approximating model at fixed bases location in the given widths. In this setting the choice of the weights is done through solving simultaneous linear algebraic equations. Considering that multiquadrics are non-local approximators, the matrices of weight systems comes nonsparse. In case the centers and widths were chosen unsuccessfully, they become ill- conditioned. Nevertheless, this paper demonstrates [14] that problem solving can be accomplished with a small set of bases due to the high accuracy of approximating functions. In these conditions the algorithm gave a lower magnitude of error and performed quicker than the finite differences method. A more general meshless approach to problem solving of heat transfer is the method recommended by A. N. Vasiliev and D. N. Tarkhov. It involves complete adjustment of all parameters of approximating solution. For this purpose, methods of unconditional nonlinear optimization are used [3, 18]. Discrete quadratic error functional is formed to achieve this. Then some set of check points is chosen in space, and the approximation model is plugged into initial equations and boundary conditions. As a rule, the error functional is smooth and in that regard it is rational to use first-order methods for accelerated search of local optimum. In a point of fact, papers [1, 2] reveal effective application of conjugate gradients non-linear method CG_DESCENT [4, 15]. An important advantage of this method is its guaranteed convergence. Alternatively, there may be considered quasi-Newton methods [5, 6], which are also very effective, but much costlier as they need more memory resources in order to keep approximated Hessian matrix (or its part, for methods with limited memory, such as L-BFGS [6]). It is necessary to understand that despite high efficiency of the optimization methods mentioned above, they are local, i.e. resulting accuracy of the solution and the time spent on its search depend heavily on initial data. Generally, the search of global optimum doesn’t bring any practical utility due to high computing costs. In a number of cases it may not even exist, particularly while considering ill-posed problems. Various approaches to global optimization, valid for parameter choice for meshless simulations, are stated in the monograph [3]. This paper presents numerical simulation of heat transfer in heterogenic multilayer structure on the macroscopic level. The paper [10] introduces an approximative model for calculating temperature balance in the layers. A finite- difference approach was used for problem solving. The computing has been exercised on implicit scheme due to physical limitation of the observed materials. In this paper we propose a hybrid method of problem solving with meshless approximation in space and finite-difference in time. 2 Problem Setting It is required to build numerical model of heat transfer in multilayer nanostructure AlAs/GaAs. It is assumed that the temperature is equal inside each layer full width, and the changes proceed in discrete steps at the interfaces. The model of the heat transfer is as follows [10]: dTi C i  i hi =  i −1 (Ti −1 − Ti ) −  i (Ti − Ti +1 ), i = 1K N , (5) dt whereC –heatcapacity, ρ –density, h –layerwidth, σ –heatconductivity, N – numberoflayers. The outer boundaries of the structure are isolated. That was achieved by zeroing of extreme values of heat conductivity:  0 =  N = 0 , then dT1 C 1  1 h1 = − i (T1 − T 2 ) , (6) dt dT N C N  N hN =  N −1 (T N −1 − T N ) . (7) dt Along with that initial temperature distribution by layers was defined: Ti ( x , 0 ) = Ti 0 ( x ) . (8) 3 Problem-Solving Method A hybrid finite-difference meshless approach in used for solving the problem. The main principle of this method involves finite-difference splitting by times and continuous approximation [1, 2] at each time layer. Thus, the algorithm performance will result in a set of continuous models, which give approximate solution at each moment of simulation. In the general case, the architecture of each approximating model may be individual, but for the simplicity, homogeneous networks with equal number of elements and identical basis function are used in this paper. The implicit finite-difference scheme serves for discretization of the first temperature derivative by time: 38 Computer Design of New Materials ___________________________________________________________________________________________ dT i k +1 ( x ) T k +1 ( x ) − T k ( x ) = . (9) dt  For dimensional approximation the radial basis functions network was used. Significantly, it is rather difficult to gain high accuracy throughout x-coordinate using RBF as infinitely smooth approximation of step functions gives a bigger error at a point of opens. In order to avoid this limitation, the bases were fixed in the centers of layers. Supposedly, the value in the center defines the temperature of the whole layer. As radial basis functions, beside of multiquadrics, the following were used for comparison:  ( x ) = exp( −  2 r 2 ) , Gaussian function, (10) 1  ( x) = ,inverse multiquadric. (11) 1 +  2r 2 Practically, normalization [13, 20] can be applied while using RBF-aproximation. In a number of cases it allows to obtain a more accurate result [1, 2, 12, 13], therefore approximative model takes the form:    ( x) . M i =1 i i u ( x) = (12)   ( x) M i =1 i The choice of parameters of the approximative model is carried out according to Kanza’s method, i.e. the centers and the widths remain fixed, and the algorithm must calculate only weight values. For this purpose, let us plug (9) in (5). After elementary transformations it becomes possible to make a tridiagonal matrix of the system of linear algebraic equations relative to the temperature of the layers T. A i , i = 1 + q i −1 + q i , A i , i −1 = − q i −1 , Ai ,i + `1 = − q i , (13)  i qi = . Ci ih In case of the classical method of finite differences the temperature calculations on the successive time layer will be conducted according the formula given below: T k +1 = A −1T k . (14) This obtained solution will be used as a standard for the accuracy assessment of meshless approximations. Numerical calculation of weight coefficients for approximating model needs the following auxiliary matrices: Bi , j = Ai ,i −1C i −1, j + Ai ,i C i , j + Ai ,i +1C i +1, j , (15) Ci, j =  i, j . In this context, the weight coefficients on the successive time step can be found in the following way: W k + 1 = B − 1T k , (16) T k +1 = CW k +1 . 4 Test Problem For testing purposes of the suggested problem solving method let us consider a heterogenic structure with N = 1000 alternating layers AlAs/GaAs, equal width for all layers h = 50 nm, thermophysical properties of materials tested: CAlAs= 424 J/(kg·К), CGaAs= 327 J/(kg·К), ρAlAs= 3730 kg/m3, ρGaAs= 5320 kg/m3 , heat conductivity on interfaces is preset as: σAlAs/GaAs= 9.6·109 Wt/(m2·К), σGaAs/AlAs= 1.92 109 Wt/(m2·К). 39 Computer Design of New Materials ___________________________________________________________________________________________ The problem solving involved normalized and unnormalized radial-basis approximators. The bases centers were 1 located in the centers of the nanostructure layers: x i = (i − 0 .5) h . The bases widths were fixed:  i = c . h Timestepwasacceptedasτ = 10-10s. Accordingly, all necessary data for setting up the system lineal algebraic equation needed for the search of weight coefficients for each time layer are given. Here are the results of the computing. The quality assessment was conducted by means of comparison with the results, obtained while using classical method of finite differences. The Figure 1, which is given below, reveals the solutions found at different moments of simulation. The numeric evaluation of accuracy needed computing of average squared deviation of the obtained solutions in the centers of layers at the considered time period. Absolute values of errors for various meshless approximations are presented in Table 1. Table 1: Absolute values of errors for various meshless approximations Approximator Error, К Multiquadric (Mq) 1.5·10-4 Gaussian function (G) 6.2·10-6 Inverse multiquadric (IMq) 3.22·10-6 Multiquadric, with normalization (Mq-N) 4.53·10-5 Gaussian function, with normalization (G-N) 3.3·10-6 Inverse multiquadric, with normalization (IMq-N) 3.52·10-6 Figure 1: Model problem solving: FDM – finite difference method; meshless methods with bases: G – Gaussian function, Мq –multiquadric, IMq – Inverse multiquadric; normalized approximations with suffix “-N It is obvious that the results shown in the diagram, which were obtained using finite differences method stiffly accurate coincide with meshless approximation. This fact demonstrates the efficiency of the method under the discussion. Table 1 shows that normalization usage allows, in some cases, to increase the accuracy of the results significantly. In fact, multiquadric has lowered the error three times, while Gaussian function has decreased its error 2 times. 5 Parallel Implementation Nowadays, one of the common ways of increasing the rate of calculation is considered to be a concurrent computation on a great number of processors. Typical usage includes multicore central processing units and graphic processing units. The algorithms, considered in this paper, utilize effectively parallelized procedures. That is the natural reason for estimating possible growth of productivity stimulated by application of concurrent computation. 40 Computer Design of New Materials ___________________________________________________________________________________________ Within the presented paper the tests were conducted in order to estimate the growth of calculation rate of the applied algorithms at parallel performance on various systems: − iMac 2012, Intel Core i5-3470S (2 cores), OSX 10.2, Clang 9; − MSI GT70, Intel Core i7-4800MQ (4 cores), Lubuntu 17.04, GCC 5.4; − IBM PowerNV, Power8 (8 cores), Ubuntu 16.04, GCC 5.4; − IBM PowerNV, Power8 (10 cores), RHEL 4.8.5, XLC_R 13.1. In order to compare productivity, synthetic tests were held, involving computation of matrix products with the help of the dgmm procedure from the standard function set for linear algebra BLAS, as execution of this function consumes a larger amount of operating hours of the program. The calculations were carried out for square matrices sized 10000х10000. Only direct calculating time was measured. For instance, time used for memory allocation was not taken into consideration. Relative results are given in Table 2. The first configuration (iMac, running single thread), was accepted as a standard. Table 2: Parallel implementation efficiency growth System (number of cores) Efficiency Growth Intel Core i5-3470S (1 core) 1 Intel Core i5-3470S (2 cores) 1.98 Intel Core i7-4800MQ (4 cores) 4.34 Power8 (8 cores) 8.94 Power8 (10 cores) 14.85 The review of the quoted results reveals that employing multicore systems allows to increase the algorithm efficiency considerably. It should be emphasized that only concurrency on CPU was under the research in this paper. The implementation of the methods mentioned above on the graphical processors will show more significant efficiency growth. 6 Conclusion The paper considered the hybrid finite difference method for heat transfer simulation in heterogenic nanostructures. The performance of the suggested algorithm was demonstrated on the example of solving a model problem. The possibility was investigated of computing rate increase by means of parallelization of calculations. It was demonstrated that usage of multiprocessing system can lead to a significant efficiency growth. The calculations were performed by Hybrid high-performance computing cluster of FRC CS RAS) [21, 22] and Shared Facility Center “Data Center of FEB RAS” (Khabarovsk) [23]. References 1. Vasilyev, A.N., Kolbin, I.S., Reviznikov, D.L.: Meshfree Computational Algorithms Based on Normalized Radial Basis Functions. Advances in Neural Networks – ISNN 2016. Lecture Notes in Computer Science, vol. 9719, pp. 583-591. (2016). DOI: 10.1007/978-3-319-40663-3 67. 2. Kolbin, I.S., Reviznikov, D.L: .Resheniye zadach matematicheskoy fiziki s ispolzovanieym normalizovannikh radialno-basisnikh setey. // «Neirokompyuteri»: razrabotka, primineniye. № 2. pp. 12-19. (2012). (In Russian) 3. Vasilyev, A.N., Tarkhov, D.A.: Neirosetevoye modelirovaniye. Principi, algorimi, prilozheniya – SPb.:Izd. Politeh. (2009). (In Russian) 4. Hager, W.W., Zhang, H.: A new conjugate gradient method with guaranteed descent and an efficient line search. SIAM Journal on Optimimization, vol. 16, pp. 170-192. (2005). DOI: 10.1137/030601880. 5. Nocedal, J., Wright, S.J.: Numerical optimization. New York: Springer Verlag, p. 634. (1999). 6. Nocedal, J.: Updating quasi-Newton matrices with limited storage. Mathematics of Computation, vol. 35, pp. 773- 782. (1980). DOI: 10.1090/S0025-5718-1980-0572855-7. 41 Computer Design of New Materials ___________________________________________________________________________________________ 7. Hardy, R.L.: Multiquadric equations of topography and other irregular surfaces. Journal of Geophysical Research, vol. 76.no. 8, pp. 1905–1915. (1971). DOI: 10.1029/JB076i008p01905. 8. Hardy, R.L.: Theory and applications of the multiquadricbiharmonic method — 20 years of discovery 1968–1988. Computers and Mathematics with Applications, vol. 19, no. 8/9, pp. 163–208. (1990). DOI: 10.1016/0898- 1221(90)90272-L. 9. Kansa, E.J.: Multiquadrics — a scattered data approximation scheme with applications to computational fluid- dynamics. 1. Surface approximations and partial derivative estimates. Computers & Mathematics with Applications, vol. 19, no. 8/9, pp. 127–45. (1990). DOI: 10.1016/0898-1221(90)90270-T. 10. Vorobyov, D.A, Khvesyuk, V.I.: Мetod rascheta nestacionarnogo nagreva nanostructur. // Nauka i obrazovaniye. – pp. 541-550. (2013). DOI: 10.7463/1013.0617255. (In Russian) 11. Chen, W.: New RBF Collocation Methods and Kernel RBF with Applications: Meshfree Methods for Partial Differential Equations. Lecture Notes in Computational Science and Engineering, vol. 26, pp. 75–86. (2002). DOI: 10.1007/978-3-642-56103-0_6. 12. Benaim, M.: On the functional approximation with normalized Gaussian units. Neural Computation, vol. 6, pp. 314-333. (1994). DOI: 10.1162/neco.1994.6.2.319. 13. Bugmann, G.: Normalized Gaussian radial basis function networks. Neurocomputing, vol. 20, no. 1/3, pp. 97-110. (1998). DOI: 10.1016/S0925-2312(98)00027-7. 14. Kansa, E.J.: Multiquadrics – a scattered data approximation scheme with applications to computational fluid dynamics-II. Solutions to hyperbolic, parabolic, and elliptic partial differential equations.Computers & Mathematics with Applications, vol. 9, no. 8/9, pp. 147-161. (1990). DOI: 10.1016/0898-1221(90)90271-K. 15. Hager, W.W., Zhang, H.: Algorithm 851: CG_DESCENT, A conjugate gradient method with guaranteed descent. ACM Transactions on Mathematical Software, no. 32, pp. 113-137. (2006). DOI: 10.1145/1132973.1132979. 16. Kansa, E.J. Sarra, S.A.: Multiquadric Radial Basis Function Approximation Methods for the Numerical Solution of Partial Differential Equations. (2009). URL: http://www.techscience.com/acm/2009/v2.html. (accessed 20.12.2016) 17. Franke, R.: Scattered data interpolation: Tests of some methods. Mathematics of Computation, vol38., no. 157, pp. 181-192. DOI: 10.1090/S0025-5718-1982-0637296-4. (1982) 18. Vasilyev, A.N., Tarkhov, D.A.: Postroyeniyepriblizhennikhneirosetvikhmodeleyporaznorodnimdannim. // Mat. Model. – 19:12. – pp. 43–51. (2007). (In Russian) 19. Haykin, S.S.: Neural Networks: A Comprehensive Foundation, 2nd Edition. Prentice-Hall, p. 936. (1999). 20. Bugmann, G., Koay, K.L., Barlow, N., Phillips M., Rodney, D.: Stable Encoding of Robot Trajectories using Normalized Radial Basis Functions: Application to an Autonomous Wheelchair. // Proc. 29th Intl. Symp. Robotics ISR'98, Birmingham, UK. – pp. 232-235. DOI: 10.1504/IJVAS.2006.01221. (1998) 21. Federal Research Center ‘Computer Science and Control’ of the Russian Academy of Sciences. Available at: http://hhpcc.frccsc.ru (accessed 09/12/2018). 22. Zatsarinny, A.A., Gorshenin, A.K., Kondrashev, V.A., Volovich, K.I., Denisov, S.A.: Toward high performance solutions as services of research digital platform. // Procedia Computer Science. Volume 150. p. 622-627. (2019) 23. Sorokin, A.A., Makogonov, S.V., Korolev, S.P.: The Information Infrastructure for Collective Scientific Work in the Far East of Russia // Scientific and Technical Information Processing. Vol. 44. Number 4. P. 302-304. (2017) 42