A simple GPU implementation of spectral-element methods for solving 3D Poisson type equations on rectangular domains and its applications

TL;DR

This paper presents a simple but extremely fast MATLAB implementation on a modern GPU, which can be easily reproduced, for solving 3D Poisson type equations using a spectral-element method and presents applications of this fast solver to solve a linear (time-independent) Schrodinger equation and a nonlinear (time-dependent) Cahn-Hilliard equation.

Abstract

It is well known since 1960s that by exploring the tensor product structure of the discrete Laplacian on Cartesian meshes, one can develop a simple direct Poisson solver with an $\mathcal O(N^{\frac{d+1}d})$ complexity in d-dimension, where N is the number of the total unknowns. The GPU acceleration of numerically solving PDEs has been explored successfully around fifteen years ago and become more and more popular in the past decade, driven by significant advancement in both hardware and software technologies, especially in the recent few years. We present in this paper a simple but extremely fast MATLAB implementation on a modern GPU, which can be easily reproduced, for solving 3D Poisson type equations using a spectral-element method. In particular, it costs less than one second on a Nvidia A100 for solving a Poisson equation with one billion degree of freedoms. We also present applications of this fast solver to solve a linear (time-independent) Schrödinger equation and a nonlinear (time-dependent) Cahn-Hilliard equation.

Figures (9)

• Figure 1: An illustration of Lagrangian $Q^k$ element and the $(k+1)\times(k+1)$ Gauss-Lobatto quadrature.
• Figure 2: Online computation time of $Q^5$ spectral-element method solving a 3D Poisson equation two hundred times. On the A100, it takes approximately only 0.8 second when solving one Poisson equation for the total number of DoFs being $1001^3$.
• Figure 3: Online computational time of $Q^5$ SEM for a Schrödinger equation with $\beta=1$ in \ref{['beta']}, solved by PCG with $(I - \Delta)^{-1}$ as the preconditioner.
• Figure 4: Comparison between $Q^5$ SEM implemented in Table \ref{['tab:3dcode']} and a second order scheme implemented by FFT in Table \ref{['code:FFT']}, for solving a Poisson equation 200 times. On A100, the FFT implementation cannot solve a problem of size $1050^3$ in MATLAB 2023, due to the larger memory cost of FFT.
• Figure 5: Comparison between $Q^5$ SEM implemented in Table \ref{['tab:3dcode']} and a second order scheme implemented by FFT in Table \ref{['code:FFT']} for solving a Schrödinger equation by PCG. On A100, the FFT implementation cannot solve a problem of size $950^3$ in MATLAB 2023, due to the larger memory cost of FFT.

Theorems & Definitions (1)

• Remark 1