ROY: Channel flow with the generalized Navier boundary condition

Purpose

• The package "Roy" simulates the 2-D incompressible channel flow with the generalized Navier boundary condition. A Legendre method was used for the spatial discretization in $y$, and a pseudo-Fourier method was used in $x$. The temporal discretization are rotational pressure correction projection schemes.

Specifications

• Name: Roy.
• Author: Feng Chen.
• Finishing date: 07/14/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{ nonhomogeneous Robin in } y, \end{aligned} \end{equation}$$ in which the generalized Navier boundary condition in $y$ is reflected: $$\begin{equation} \begin{aligned} &\boldsymbol{v} \cdot \boldsymbol{n} = 0, && \text{on } \partial\Omega, \\ & \big[ \boldsymbol{\sigma} \boldsymbol{n}+ l(\phi) \boldsymbol{v}_{\text{s}} \big]\times \boldsymbol{n}=g(x,y,t), && \text{on } \partial\Omega. \end{aligned} \end{equation}$$
• Parameters: $$\begin{equation} \Omega =(-\pi,\pi)\times (-1,1), \, T=0.2, \, N_x = 64, \, N_y = 32 . \end{equation}$$
• Exact solution and input functions: $$\begin{equation} \begin{aligned} & v_1 = \pi \sin(2\pi y)\sin^2(x) \sin(t), \\ & v_2 = - \sin^2(\pi y)\sin(2x) \sin(t), \\ & p = \cos(x) \cos(\pi y) \sin(t), \end{aligned} \end{equation}$$ such that $$\begin{equation} \bar{l} v_1 \pm \eta v_1' = \pm 2\pi^2 \eta \sin^2(x) \sin(t), \quad y = \pm 1, \end{equation}$$ where $\bar{l}=4$ and $\eta=3$. And $\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:

• 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 rotational pressure 
!! correction scheme.
!<
! external force.
call Roy_Utility_Time_F(RF, is*RoyT%dt)

! convection of v.
call IonPhotonSpace_Transform_D2_Many(2, RoySI, RoySP, &
 RF%vp(:,:,:,co-1), RF%vs, IonPhotonP2S)
call IonPhotonSpace_Convection_D2(RoySI, RoySP, RF%vs, &
 RF%conv(:,:,:,co-1), RF%gp, RF%pfi)

! gradient of p.
call IonPhotonSpace_Transform_D2_Many(1, RoySI, RoySP, &
 RF%pp(:,:,co-1), RF%ps, IonPhotonP2S)
call IonPhotonSpace_Gradient_D2(RoySI, RoySP, RF%ps, RF%gp)
call IonPhotonSpace_Transform_D2_Many(2, RoySI, RoySP, &
 RF%gp, RF%gp, IonPhotonS2P)

! cp, cm, and fd for v.
call Roy_Utility_Slip(RF, co, is*RoyT%dt)
call IonSpace_Transform_D1_Many(1, RoySI, RF%cp, RF%cp, IonP2S)
call IonSpace_Transform_D1_Many(1, RoySI, RF%cm, RF%cm, IonP2S)
call IonQuarkSpace_RHSData(RoySI, RoySQv1(co), &
 RF%cp, RF%cm, RF%fd)

! solve for v.
RF%vp(:,:,:,co) = 0d0
do io = 0, co - 1
 a = RoyT%a(co,io)
 b = RoyT%b(co,io)
 RF%vp(:,:,:,co) = RF%vp(:,:,:,co) - (a/RoyT%dt) * &
  RF%vp(:,:,:,io) - b * RF%conv(:,:,:,io)
end do
RF%vp(:,:,:,co) = RoyRe*RF%vp(:,:,:,co) - RF%gp + &
  RoyB*RF%pfi + RF%f
call IonPhotonSpace_Transform_D2_Many(2, RoySI, RoySP, &
 RF%vp(:,:,:,co), RF%vs, IonPhotonP2S)
call IonQuarkSpace_Solve(RoySI, RoySQv1(co), RF%fd, RF%vs(:,:,1))
call IonAtomSpace_Solve_D2(1, RoySI, RoySAv2, RoyCaV2(co), &
 RoyCb1V2(co), RoyCb2V2(co), RF%vs(:,:,2))
call IonPhotonSpace_Transform_D2_Many(2, RoySI, RoySP, RF%vs, &
 RF%vp(:,:,:,co), IonPhotonS2P)

! divergence of v.
call IonPhotonSpace_Divergence_D2(RoySI, RoySP, RF%vs, &
 RF%div, RF%gp)

! update v and p.
call IonPhotonSpace_Gradient_D2(RoySI, RoySP, RF%ps, RF%gp)
call IonPhotonSpace_Transform_D2_Many(2, RoySI, RoySP, RF%gp, &
 RF%gp, IonPhotonS2P)
RF%vp(:,:,:,co) = RF%vp(:,:,:,co) - (RoyT%dt/ &
 (RoyRe*RoyT%a(co,co)) ) * RF%gp

call IonPhotonSpace_Transform_D2_Many(1, RoySI, RoySP, &
 RF%div, RF%div, IonPhotonS2P)
RF%pp(:,:,co) = RF%pp(:,:,co-1) + RF%pp(:,:,co) - &
 RoyEta * RF%div