GOLD: Spectral methods for 3-D driven cavity flows

Purpose

• The package "Gold" simulates the 3-D incompressible driven cavity flow. The Legendre-Galerkin method was used for the spatial discretization, and the first- and second-order rotational pressure correction projection methods were used for the temporal discretization.

Specifications

• Name: Gold.
• Author: Feng Chen.
• Finishing date: 07/06/2011.
• Languages: Fortran 90.
• Required libraries: BLAS, LAPACK.

Simple Example

• Equation: \begin{equation} \begin{aligned} & \boldsymbol{v}_t + \boldsymbol{v} \cdot \nabla \boldsymbol{v} = \nu \Delta \boldsymbol{v} - \nabla p + \boldsymbol{f}, \\ & \nabla \cdot \boldsymbol{v} = 0, \\ & \boldsymbol{v}(x,y,z,0) = \boldsymbol{v}_0(x,y,z), \\ & \boldsymbol{v}|_{\partial \Omega} = \boldsymbol{0}. \end{aligned} \end{equation}
• Parameters: \begin{equation} \Omega =(-1,1)^3, \, T=1, \, \delta t = 0.01, \, N_x = N_y= N_z=64, \, \nu = 1. \end{equation}
• Exact solution and input functions: \begin{equation} \begin{aligned} & v_1(x,y,z,t) = 2\pi \sin^2 (\pi x) \sin(2\pi y) \sin(2\pi z) \sin(t), \\ &v_2(x,y,z,t) = -\pi \sin(2\pi x) \sin^2 (\pi y) \sin(2\pi z) \sin(t),\\ &v_3(x,y,z,t) = -\pi \sin(2\pi x) \sin (2\pi y) \sin^2(\pi z) \sin(t),\\ & p(x,y,z,t) = \cos(\pi x) \cos(\pi y) \cos (\pi z) \sin(t), \end{aligned} \end{equation} then $\boldsymbol{f}(x,y,t)$ is calculated accordingly.

Quick Start

• Compiling and running:

  cd ./Gold
  make library
  gfortran Gold_Main.cu -llibrary -llapack -lblas
  ./a.out
  

• Output:

   
  Second-order scheme was used. 
   At step           20 Error of v =  1.71716801470710173E-004
   At step           20 Error of p =  4.30641138903224419E-002
   At step           40 Error of v =  1.61822030845909292E-004
   At step           40 Error of p =  7.41918977991906092E-002
   At step           60 Error of v =  1.45463465091464418E-004
   At step           60 Error of p =  0.10235141260382691     
   At step           80 Error of v =  1.23321771533398253E-004
   At step           80 Error of p =  0.12646838834271690     
   At step          100 Error of v =  9.62650234262930784E-005
   At step          100 Error of p =  0.14556457932094957     
  ...
  

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

References

• Guermond, J. L. and Shen, J. On the error estimates for the rotational pressure-correction projection methods, Mathematics of Computation, No. 248 (2003).
• Shen, J. Efficient spectral-Galerkin method I. Direct solvers of second-and fourth-order equations using Legendre polynomials, SIAM Journal on Scientific Computing, Vol. 15, No. 6 (1994).

Code Highlights

!<
!! A typical procedure in the rotational pressure correction
!! projection scheme
!<

!! calculate the nonlinear term and the gradient of the pressure.   
call GoldSpace_Calculate_Convection(SV, B%velocity(:,:,:,:,co-1)+ &
 GoldVic, B%conv(:,:,:,:,co-1))
call GoldSpace_Calculate_Gradient(SP, B%pressure(:,:,:,co-1), B%gradp) 

!! form the right side of the momentum equation.   
B%velocity(:,:,:,:,co) = 0d0
do io = 0, co - 1
 B%velocity(:,:,:,:,co) = B%velocity(:,:,:,:,co) - &
  (T%a(co,io)/T%dt) * B%velocity(:,:,:,:,io) - &
  T%b(co,io) * B%conv(:,:,:,:,io)
end do
B%gradp2 = 0d0  
B%gradp2(0:SP%Np(1),0:SP%Np(2),0:SP%Np(3),1:3) = B%gradp
B%velocity(:,:,:,:,co) = B%velocity(:,:,:,:,co) - B%gradp2 + &
 B%exfr + GoldVisco*GoldVicdd   

!! solve for the velocity field. 
call GoldSpaceTime_Solve_velocity(co, SV, B%velocity(:,:,:,:,co))  

!! calculate divergence of the velocity field. 
call GoldSpace_Calculate_Divergence(SV, &
 B%velocity(:,:,:,:,co), B%divv) 

!! form the right side of the projection equation.  
B%divv2 = B%divv(0:SP%Np(1),0:SP%Np(2),0:SP%Np(3))
B%pressure(:,:,:,co) = -(T%a(co,co)/T%dt) * B%divv2    

!! solve for the (auxiliary) B%pressure. 
call GoldSpaceTime_Solve_pressure(co, SP, B%pressure(:,:,:,co))   

!! update the velocity and the pressure. 
call GoldSpace_Calculate_Gradient(SP, B%pressure(:,:,:,co), B%gradp)
B%gradp2 = 0d0 
B%gradp2(0:SP%Np(1),0:SP%Np(2),0:SP%Np(3),1:3) = B%gradp
B%velocity(:,:,:,:,co) = B%velocity(:,:,:,:,co) - &
 (T%dt/T%a(co,co)) * B%gradp2
B%pressure(:,:,:,co) = B%pressure(:,:,:,co-1) + &
 B%pressure(:,:,:,co) - GoldVisco * B%divv2