SILVER: Fourier-spectral-element for channel flow

Purpose

• The package "Silver" simulates the 2-D incompressible channel flow. A Fourier psuedospectral method was used for the spatial discretization in $x$, and a modal spectral element method was used for $y$. The temporal discretization are rotational pressure correction projection schemes. It is partially based on packages "Sea" and "Ion".

Specifications

• Name: Silver.
• Author: Feng Chen.
• Finishing date: 07/12/2012.
• 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,0) = \boldsymbol{v}_0(x,y), \\ & \boldsymbol{v} \text{ periodic in } x, \\ & \boldsymbol{v} \text{ homogeneous Dirichlet in } y. \end{aligned} \end{equation}
• Parameters: \begin{equation} \Omega =(-\pi,\pi)\times (-1,1), \, T=10, \, \nu = 1,\, \delta t = 0.001, \, N_x = 64, \, N_y = 10, \, K_y = 10. \end{equation}
• Exact solution and input functions: \begin{equation} \begin{aligned} & v_1(x,y,t) = \pi \sin(2\pi y) \sin^2 (x) \sin(t), \\ &v_2(x,y,t) = \sin(2 x) \sin^2 (\pi y) \sin(t),\\ & p(x,y,t) = \cos(x) \cos(\pi y) \sin(t), \end{aligned} \end{equation} then $\boldsymbol{f}(x,y,t)$ is calculated accordingly.

Quick Start

• Compiling and running:

  cd ./Silver
  make library
  gfortran Silver_Main.cu -llibrary -llapack -lblas
  ./a.out
  

• Output:

  Order 1 scheme was used.
   Max error of v at step           0 :   0.0000000000000000     
   Max error of p at step           0 :   0.0000000000000000     
   Max error of v at step        1000 :  1.36078134067213474E-004
   Max error of p at step        1000 :  3.37562858666562438E-003
   Max error of v at step        2000 :  9.62240734776464990E-005
   Max error of p at step        2000 :  2.67380732043620561E-003
   Max error of v at step        3000 :  4.79183164863519750E-005
   Max error of p at step        3000 :  1.04866367561439605E-003
   Max error of v at step        4000 :  1.44443446522313224E-004
   Max error of p at step        4000 :  3.72881892945142734E-003
   Max error of v at step        5000 :  6.84469621439198761E-005
   Max error of p at step        5000 :  1.98585017908048389E-003
   Max error of v at step        6000 :  8.26089424067338873E-005
   Max error of p at step        6000 :  2.00780300921105104E-003
   Max error of v at step        7000 :  1.39257211530317837E-004
   Max error of p at step        7000 :  3.62196919613311241E-003
   Max error of v at step        8000 :  4.13692986276359420E-005
   Max error of p at step        8000 :  1.07504147050024867E-003
   Max error of v at step        9000 :  1.09539686892368149E-004
   Max error of p at step        9000 :  2.76938279571875556E-003
   Max error of v at step       10000 :  1.27007787214084011E-004
   Max error of p at step       10000 :  3.47656709161936472E-003
  ...
  

• 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).

Code Highlights

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

!! calculate the nonlinear convection term 
!! and the gradient of SS%p. 
call SilverSpace_Calculate_Convection(SI, SV, SS%v(:,:,:,co-1)+&
 SilverVic, SS%cv(:,:,:,co-1), SilverWk2, SilverWk3)
call SilverSpace_Calculate_Gradient(SI, SP, SS%p(:,:,co-1), SS%gp)

!! form the right side of the v equation.    
SS%v(:,:,:,co) = 0d0
do io = 0, co - 1
 SS%v(:,:,:,co) = SS%v(:,:,:,co) - &
  (T%a(co,io)/T%dt) * SS%v(:,:,:,io) - &
   T%b(co,io) * SS%cv(:,:,:,io)  
end do
SS%v(:,:,:,co) = SS%v(:,:,:,co) + SilverLambda*SS%fpi(:,:,:) 
SS%v(:,:,:,co) = SS%v(:,:,:,co) - SS%gp + &
 SS%exfr + SilverVisco*SilverVicdd  

!! solve for v. 
call SilverSpaceTime_Solve_Velocity(co, SI, SV, SS%v(:,:,:,co))    
 
!! calculate divergence of v. 
call SilverSpace_Calculate_Divergence(SI, SV, SS%v(:,:,:,co), SS%dv)

!! form the right side of the projection equation.    
SS%p(:,:,co) = -(T%a(co,co)/T%dt) * SS%dv 

!! solve for the (auxiliary) SS%p. 
call SilverSpaceTime_Solve_Pressure(co, SI, SP, SS%p(:,:,co))     

!! update the SS%v and SS%p. 
call SilverSpace_Calculate_Gradient(SI, SP, SS%p(:,:,co), SS%gp) 
SS%v(:,:,:,co) = SS%v(:,:,:,co) - (T%dt/T%a(co,co)) * SS%gp 
SS%p(:,:,co) = SS%p(:,:,co-1) + SS%p(:,:,co) - SilverVisco * SS%dv     

!! calculate the nonhomogeneous velocity vpt. 
SilverWk2 = SS%v(:,:,:,T%ot) + SilverVic
do i = 1, 2
 call SilverSpace_Transform(SI, SV, SilverWk2(:,:,i), &
  SS%vpt(:,:,i), SilverS2P)    
end do