LION: Scalable Legendre-Galerkin methods for elliptic systems

Purpose

• The package "Lion" provides a Legendre-Galerkin method for solving systems of coupled elliptic equations with the Dirichlet, Neumann, or Robin boundary conditions on a 2-D or 3-D rectangular domain.

Specifications

• Name: Lion.
• Author: Feng Chen.
• Finishing date: 04/08/2012.
• Language: Fortran 90.
• Required libraries: BLAS, LAPACK.

Simple Example

• Equation: $$\begin{equation} \begin{aligned} &\alpha u + \beta v - \Delta u = f, &&\quad \text{in } \Omega, \\ &\tilde \alpha u + \tilde \beta v - \Delta v = g, &&\quad \text{in } \Omega, \\ &u = v = 0, && \quad \text{on } \partial \Omega. \end{aligned} \end{equation}$$
• Parameters: $$\begin{equation} \Omega = (-1,1)^2 , \quad \alpha =2, \quad \beta=1, \quad \tilde \alpha = 1, \quad \tilde \beta=2. \end{equation}$$
• Exact solution and input functions: $$\begin{equation} \begin{aligned} & u(x,y) = \sin(\pi x) \sin(2\pi y), \\ & v(x,y) =\sin(2\pi x) \sin(3\pi y). \end{aligned} \end{equation}$$ then $f(x,y)$ and $g(x,y)$ are calculated accordingly.

Quick Start

• Compiling and running:

  cd ./Lion
  make library
  gfortran tst_Lion.f90 -llibrary -llapack -lblas
  ./a.out
  

• Output:

   
  Nx = 256, Ny = 128
  computing time = 0.10900000 sec
  max error of u = 0.34E-14
  max error of v = 0.37E-14
  ...
  

• CPU: Intel(R) Xeon(R) CPU X5550 @2.67GHz.
• OS: CentOS release 6.4 (Final).
• Compiler: gfortran 4.5.1.

References

• Feng Chen and Jie Shen. Efficient spectral-Galerkin methods for systems of coupled second-order equations and their applications, Journal of Computational Physics, Volume 231, Issue 15, 5016-5028, (2012).

Code Highlights

!! L:      number of equations in the system
!! n:      number of quadrature points
!! q:      number of basis
!! a,b,d:  contruction of basis
!!  ca,cb:  equation coefficients
!! fs:     solution vector of size nnL
!! E:      orthnormal basis
!! lambda: eigenvalues
!! s:      stiffness matrix

!! calculate the right hand side (inner products)
do m = 1, L 
 call D2_Legen_Gal_Opti_RHS(nx, ny, qx, qy, &
  ax, ay, bx, by, dx, dy, fs(:,:,m))
end do

!! transform to the orthogonal space (diagonalization)
do m = 1, L
 fs(0:qx,0:qy,m) = matmul(transpose(Ex), fs(0:qx,0:qy,m))
 fs(0:qx,0:qy,m) = matmul(fs(0:qx,0:qy,m), Ey)
end do

!! inversion in the spectral space (pointwise inversion)
do i = 0, qx
 do j = 0, qy
  !! wave number
  cw = Lambdax(i)*Lambday(j)*ca  &
   + cb1*sx(i)*Lambday(j)  &
   + cb2*Lambdax(i)*sy(j)      
  !! for the Neumann boundary condition
  if (maxval(abs(cw)) .eq. 0d0) then
   do m = 1, L
    cw(m,m) = 1d0
   end do      
  end if
  !! inversion
  call DGESV( L, 1, cw, L, IPIV, fs(i,j,:), L, INFO )
 end do
end do

!! transform to the solution space (substitution)
do m = 1, L
 fs(0:qx,0:qy,m) = matmul(Ex, fs(0:qx,0:qy,m))
 fs(0:qx,0:qy,m) = matmul(fs(0:qx,0:qy,m), transpose(Ey))
end do

!! from basis form to spectral form
do m = 1, L 
 call D2_Legen_Gal_Opti_BasisTransform(qx, qy, nx, ny, &
  ax, ay, bx, by, dx, dy, fs(:,:,m))
end do