Numerical solution of elastodynamics equations using finite differences and dynamic calls on GPUs

Mathematical models based on elastodynamics equations are often used to study earthquake dynam- ics [1], [2]. These models lead to the production of strikingly non-trivial patterns. Therefore, their numerical solution often requires high spatial resolution to capture the detailed parameters related to the phenomena, and as a consequence, long computation times are often required when using a se- rial implementation of a numerical scheme. Parallel computation can dramatically improve the time and efficiency of some numerical methods such as finite difference algorithms, which are relatively simple to implement and apply to different models. For applied scientists involved in setting up re- alistic experiments, the possibility of running fast comparative simulations using simple algorithms implemented on affordable processors is of primary interest, and that is where Graphical Processing Units (GPUs) can excel. Parallel computing based on modern GPUs has the advantage of a high performance at relatively low energy and monetary costs. The codes used to study the performance of GPUs presented in this article were programmed us- ing CUDA FORTRAN [3], [4]. We use three programming techniques (standard CUDA, dynamic calls and streams) to solve the elastodynamics equations, for which the numerical solution is obtained using a second order finite difference method for the spatial and time discretizations, and we com- pared the speed up and efficiency related to each code. The choice of time-explicit algorithms is due to their greater ease of implementation and performance, and despite their limitations related to reduced stability properties. The main target of this work it is a preliminary study based on the bino- mial finite-difference method and GPU computing that will be the base to a more complex numerical scheme that will merge Discontonuos Galerkin method for the treatment of irregular boundary and finite difference method for computing on the bulk of the treated domains