AMBER: GPU-accelerated 2-D Chebyshev collocation methods for elliptic equations

Purpose

• The package "Amber" provides a GPU-accelerated Chebyshev collocation method for solving the single elliptic equation with Dirichlet, Neumann, or Robin boundary conditions on a 2-D rectangular domain.

Specifications

• Name: Amber.
• Author: Feng Chen.
• Finishing date: 04/15/2012.
• Languages: CUDA C++, Fortran 90.
• Required libraries: BLAS, LAPACK, CUBLAS.

Simple Example

• Equation: \begin{equation} \begin{aligned} & \alpha u -(\beta_1 u_{xx} + \beta_2 u_{yy}) = f, &&\quad \text{in } \Omega, \\ & au + b\frac{\partial u}{\partial \boldsymbol{n}} = c, && \quad \text{on } \partial \Omega. \end{aligned} \end{equation}
• Parameters: \begin{equation} \Omega = (-1,1)^2 , \quad \alpha =2, \quad \beta_1=3, \quad \beta_2=4, \quad a = 1, \quad b=1. \end{equation}
• Exact solution and input functions: \begin{equation} u(x,y) = e^{x+y}, \end{equation} then $f(x,y)$ and $c(x,y)$ are calculated accordingly.

Quick Start

• Compiling and running:

  cd ./Amber
  make library
  nvcc Amber_Main.cu -llibrary -llapack -lblas
  ./a.out
  

• Output:

  N = 2^3, GPU time = 1.5e-04
  N = 2^4, GPU time = 4.7e-05
  N = 2^5, GPU time = 5.3e-05
  N = 2^6, GPU time = 6.8e-05
  ...
  

• CPU: Intel(R) Xeon(R) CPU X5630 @2.53GHz.
• GPU: Nvidia Tesla M2050.
• OS: CentOS release 6.4 (Final).
• Compiler: nvcc 4.2.9, gcc/gfortran 4.5.1.

References

• Feng Chen and Jie Shen. A GPU parallelized spectral method for elliptic equations in rectangular domains, Journal of Computational Physics, Volume 250, 555-564, (2013).

Code Highlights

// 
// Point-wise inversion in the frequency space
// Input:  alpha, betax, betay (coefficients)
//         lambdax, lambday (eigenvalues)
//         nx, ny (number of quadrature points)
//         f (input function)
// Output: f (output function)
//
__global__ void Mercury_Kernel(double alpha, 
 int nx, double betax, double* lambdax, 
 int ny, double betay, double* lambday, 
 double* f) {
 
 // index of the current thread
 int i = threadIdx.x + blockIdx.x * blockDim.x;
 int j = threadIdx.y + blockIdx.y * blockDim.y;
 
 // index of the element to be inversed
 int i2;
 
 // wave number
 double sigma;
 
 if ( (i < nx) && (j < ny) ) {  
  sigma = alpha - betax*lambdax[i] - betay*lambday[j];  
  // for homogeneous Neumann condition and alpha=0 
  if (fabs(sigma) < 1e-8) {
   sigma = 1.0;
  }  
  // point-wise inversion
  i2 = c2(i,j,nx);
  f[i2] = f[i2]/sigma;
 }        
}