A GPU parallelized spectral method for elliptic equations in rectangular domains

Abstract

We design and implement a polynomial-based spectral method on graphic processing units (GPUs). The key to success lies in the seamless integration of the matrix diagonalization technique and the new generation CUDA tools. The method is applicable to elliptic equations in rectangular domains with general boundary condition. We show remarkable speedups of up to 15 times in the 2-D case and more than 35 times in the 3-D case.

Introduction

General-purpose computing on graphics processing units (GPGPU) has drawn much attention recently from the scientific computing community. The new list of top 500 supercomputers shows that more than 10% of them are now powered by NVIDIA Tesla GPUs (cf. [4]). Thanks to the high core density and the wide vector width SIMD architecture, using GPGPU can yield performance that is very hard for a conventional CPU to achieve, especially for high and regular throughput workloads (cf. [3]). Therefore, researchers in computational sciences are more and more interested in exploring efficient parallel strategies on GPUs.

There exist many successful GPU implementations of numerical methods for partial differential equations, to name a few, fast multipole methods [10], nodal discontinuous Galerkin methods [12], finite difference methods [6], finite element methods [24], [14], [13], Fourier spectral methods [9], [5], etc. However, to the best of the authors’ knowledge, not much effort, if not at all, are devoted to implementing, on GPUs, non-Fourier based spectral methods for problems with non-periodic boundary conditions. The aim of this paper is to exploit how to efficiently implement spectral methods with GPGPU. In particular, we shall design and implement an efficient matrix diagonalization method with spectral discretizations on the GPU, for solving elliptic and parabolic type equations with non-periodic boundary conditions on a d-dimensional rectangular domain ($d=2\text{,}3$).

It is well-known that, for separable elliptic equations, one can use the so-called matrix diagonalization technique (cf. [16], [11], [18]) to solve the linear systems from a spectral discretization in O($N^{d+1}$) operations, where N is the number of points in each dimension. For time-dependent problems, these solvers are called repeatedly. Fortunately, the core part of this matrix diagonalization process is matrix–matrix multiplications, which can be efficiently parallelized on the GPU. More precisely, matrix–matrix multiplications can be efficiently addressed by the new generation CUDA tool, CUBLAS.

We found in the numerical experiments that the speedup of using GPU vs. CPU goes up as N increases. For example, the speedup reaches a factor of 15 with $N=2^{12}$ in the 2-D case, and a factor of 37 with $N=2^9$ in the 3-D case, where N is the number of points in each direction. In addition, computing times of different parts scale consistently with the complexity. In our algorithm, the throughput is regular and the global memory access is completely coalesced, indicating a good performance on the GPU.

To initialize the matrix diagonalization procedure, we need to first compute the eigenvalue and eigenvectors of 1-D problems on the host, then ship these eigenvalues and eigenvectors to the device because they will be used repeatedly in a typical time marching scheme for evolutionary problems.

This paper is organized as follows. Details of the matrix diagonalization method are described in Section 2 using the spectral-collocation method as an example. They provide a background for an in-depth discussion in Section 3 for the GPU parallelization of the same algorithm. In Section 4, we present some benchmark results for solving model elliptic equation in both 2-D and 3-D, and for solving 2-D multiphase flows based on the coupled system of Navier–Stokes equation and Allen–Cahn equation.