Porting a high-order finite-element earthquake modeling application to NVIDIA graphics cards using CUDA

Abstract

We port a high-order finite-element application that performs the numerical simulation of seismic wave propagation resulting from earthquakes in the Earth on NVIDIA GeForce 8800 GTX and GTX 280 graphics cards using CUDA. This application runs in single precision and is therefore a good candidate for implementation on current GPU hardware, which either does not support double precision or supports it but at the cost of reduced performance. We discuss and compare two implementations of the code: one that has maximum efficiency but is limited to the memory size of the card, and one that can handle larger problems but that is less efficient. We use a coloring scheme to handle efficiently summation operations over nodes on a topology with variable valence. We perform several numerical tests and performance measurements and show that in the best case we obtain a speedup of 25.

Introduction

Over the past several years, non-graphical applications ported to the GPUs have steadily grown more numerous [25]. Many applications have been adapted to the GPU using one of several specialty languages. These range from graphical languages such as Cg, HLSL, GLSL, to languages that abstract away the graphics component to ease the implementation of non-graphical code, such as Brook, Scout [20], Sh, Rapidmind, and CUDA. CUDA offers significant innovation and gives users a significant amount of control over a powerful streaming/SIMD computer. Since CUDA became available in 2007, the number of General-purpose Processing on Graphical Processing Units (GPGPU) applications has drastically increased. CUDA is being used to study physical problems as diverse as molecular dynamics [2], fluid dynamics [4], or astrophysical simulations [24].

Our research group has developed a software package, called SPECFEM3D_GLOBE, for the three-dimensional numerical simulation of seismic wave propagation resulting from earthquakes in the Earth based on the so-called spectral-element method (SEM) [16], [19], [17]. The SEM is similar to a finite-element method with high-degree polynomial basis functions. The mesh of elements required to cover the Earth is usually so large that calculations are very expensive in terms of CPU time and require a large amount of memory, at least a few terabytes. We therefore usually turn to parallel programming with MPI [17], [15], OpenMP, or both. Typical runs require a few hundred processors and a few hours of elapsed wall-clock time. Current large runs require a few thousand processors, typically 2000–4000, and two to five days of elapsed wall-clock time [17], [15]. It is therefore of interest to try to speed up the calculations on each node of such a parallel machine by turning to GPGPU. We will see in the next section that a significant amount of work has been devoted in the community to porting low-order finite-element codes to GPU, but to our knowledge we are the first to address the issue of porting a more complex higher-order spectral-like finite-element technique.

Rather than attempt to solve the fully resolved problem using both MPI and the GPU simultaneously, which we will consider in a future study, we choose to investigate the potential for speedup on a single CPU, with acceleration provided solely by the GPU. Such an implementation can be very useful in the absence of a GPU cluster for quick parameter studies, parametric explorations, and code development. In the future the experience gained by this exercise will provide the proper intuition to accelerate an MPI version of the code using the graphics processors. Our seismic wave propagation application is written in Fortran95 + MPI, but we wrote a serial version in C for the tests to facilitate interfacing with CUDA, which is currently easier from C.

A major issue in all finite-element codes, either low or high order, is that dependencies arise from the summation of elemental contributions at those global points shared among multiple mesh elements, requiring atomic operations. We will show that use of a mesh coloring technique allows us to overcome this difficulty and obtain a speedup of 25x, consistent with speedups obtained by many other researchers for other realistic applications in different fields.

The remainder of the paper is organized as follows: We begin with a description of similar problems ported to the GPU in Section 2. The serial algorithm is discussed in Section 3, followed by the CUDA implementation of the spectral-element algorithm in Section 4. Within this section, we discuss an implementation that runs fully on the GPU, including the coloring scheme that is needed to avoid atomic operations. We also describe a version of our algorithm capable of solving larger problems that cannot fit fully on the GPU. Optimizations we considered to improve efficiency are treated in Section 5. We provide some numerical code validation in Section 6, and analyze some performance results in Section 7. We conclude in Section 8.