IRON: Spectral methods for 2-D driven cavity flows

Purpose

• The package "Iron" simulates the 2-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: Iron.
• Author: Feng Chen.
• Finishing date: 06/30/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,0) = \boldsymbol{v}_0(x,y), \\ & \boldsymbol{v}|_{\partial \Omega} = \boldsymbol{0}. \end{aligned} \end{equation}$$
• Parameters: $$\begin{equation} \Omega =(-1,1)^2, \, T=1, \, \delta t = 0.001, \, N_x = N_y= 128, \, \nu = 1. \end{equation}$$
• Exact solution and input functions: $$\begin{equation} \begin{aligned} & v_1(x,y,t) = \pi \sin(2\pi y) \sin^2 (\pi x) \sin(t), \\ &v_2(x,y,t) = \pi \sin(2\pi x) \sin^2 (\pi y) \sin(t),\\ & p(x,y,t) = \cos(\pi 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 ./Iron
  make library
  gfortran Iron_Main.cu -llibrary -llapack -lblas
  ./a.out
  

• Output:

   
  Second-order scheme was used.
  ---------            0        init -----------------
   Error of v =   0.0000000000000000     
   Error of p =   0.0000000000000000     
   ---------           1         250 -----------------
   Error of v =  2.00242220500198424E-006
   Error of p =  3.29543665066828195E-003
   ---------           2         500 -----------------
   Error of v =  1.81394643174368083E-006
   Error of p =  6.30056824367714041E-003
   ---------           3         750 -----------------
   Error of v =  1.51268851000740494E-006
   Error of p =  8.91364994804177968E-003
   ---------           4        1000 -----------------
   Error of v =  1.11737921600563328E-006
   Error of p =  1.09726145825772559E-002
  ...
  

• 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 convection term 
!! and the gradient of pressure. 
call IronSpace_Calculate_Convection(SV, B%velocity(:,:,:,co-1)+&
 IronVic, B%conv(:,:,:,co-1))
call IronSpace_Calculate_Gradient(SP, B%pressure(:,:,co-1), B%gradp) 

!! test grad(p), can be commented. 
do i = 1, 2
 call IronSpace_Transform(SP, B%gradp(:,:,i), &
  IronGradPphy(:,:,i), IronS2P)  
end do    

!! form the right side of the velocity 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) + &
  IronLambda*T%b(co,io)*B%pfi(:,:,:,io)
end do 
B%gradp2 = 0d0  
B%gradp2(0:SP%Np(1),0:SP%Np(2),1:2) = B%gradp
B%velocity(:,:,:,co) = B%velocity(:,:,:,co) - B%gradp2 + &
 B%exfr + IronVisco*IronVicdd   

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

!! calculate divergence of the velocity field. 
call IronSpace_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))
B%pressure(:,:,co) = -(T%a(co,co)/T%dt) * B%divv2    

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

!! update the velocity and pressure. 
call IronSpace_Calculate_Gradient(SP, B%pressure(:,:,co), B%gradp)
B%gradp2 = 0d0
B%gradp2(0:SP%Np(1),0:SP%Np(2),1:2) = 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) - IronVisco * B%divv2