erturk_2.f

C****************************************************************************C
C                                                                            C
C     This program solve the Navier-Stokes (NS) equations using the          C
C     method proposed by Erturk et al. ["Numerical Solutions of 2-D          C
C     Steady Incompressible Driven Cavity Flow at High Reynolds Numbers",    C
C     International Journal for Numerical Methods in Fluids (2005),          C
C     Vol 48, pp 747-774]. The solution is spatially second order            C
C     accurate. Please refer to this paper for details.                      C
C     This code is written by Prof. Ercan Erturk.                            C
C     Visit http://www.cavityflow.com                                        C
C                                                                            C
C********************************************C*******************************C
C                                            C
C     s(i,j) ==> streamfunction variable     C
C     v(i,j) ==> vorticity variable          C
C     x(i)   ==> x-coordinate                C
C     y(j)   ==> y-coordinate                C
C     dh     ==> grid spacing                C
C     Re     ==> Renolds Number              C
C                                            C
C********************************************C
      

      program main
      
      implicit double precision (a-h,o-z)
      
      parameter(N=128)
      
      common / flow variables /
     > s(0:N,0:N),v(0:N,0:N),
     > s_old(0:N,0:N),v_old(0:N,0:N)
      common / geometry /
     > x(0:N),y(0:N)
      common / other /
     > asx(0:N,0:N),bsx(0:N,0:N),csx(0:N,0:N),
     > asy(0:N,0:N),bsy(0:N,0:N),csy(0:N,0:N),
     > avx(0:N,0:N),bvx(0:N,0:N),cvx(0:N,0:N),
     > avy(0:N,0:N),bvy(0:N,0:N),cvy(0:N,0:N),
     > s_rhs_1(0:N,0:N),v_rhs_1(0:N,0:N),
     > s_rhs_2(0:N,0:N),v_rhs_2(0:N,0:N),
     > ds(0:N,0:N),dv(0:N,0:N),
     > ff(0:N,0:N),gg(0:N,0:N)


       Re=1000.d0

       dh=1.0d0/dble(N)

      do 1 k=0,N
       x(k)=dble(k)/dble(N)
       y(k)=dble(k)/dble(N)
 1    continue

C     Time step
       alpha=1.6d0
       dt=alpha*dh*dh

C     Initial guess
c     Note that when homogeneous initial guess is used, in the first iteration 
c     residual_3 is indeterminate. In order to avoid this, instead of using 
c     exactly zero initial values, at interior points use very very small 
c     numbers that could be considered as zero.
      do 2 k=0,N
      do 22 j=0,N
       s(k,j)=1.d-32
       v(k,j)=1.d-32
 22   continue
       s(0,k)=0.d0
       v(0,k)=0.d0
       s(N,k)=0.d0
       v(N,k)=0.d0
       s(k,0)=0.d0
       v(k,0)=0.d0
       s(k,N)=0.d0
       v(k,N)=0.d0
 2    continue

C     Record the CPU time at start
       call cpu_time(time_start)

C     Start the iterations
      do 999 iteration=1,1000000

C     Update old variables
      do 5 i=1,N-1
      do 5 j=1,N-1
       s_old(i,j)=s(i,j)
       v_old(i,j)=v(i,j)
 5    continue

C     SOLVE THE STREAMFUNCTION EQUATION        

c     Construction of matrices for streamfunction equation
      do 100 i=1,N-1
      do 100 j=1,N-1
       asx(i,j)=-dt/dh**2.
       bsx(i,j)=1.d0+dt*2.d0/dh**2.
       csx(i,j)=-dt/dh**2.

       asy(i,j)=-dt/dh**2.
       bsy(i,j)=1.d0+dt*2.d0/dh**2.
       csy(i,j)=-dt/dh**2.
 100  continue

c     Calculate the added terms to RHS
      do 101 j=1,N-1
       s_rhs_1(0,j)=0.d0
      do 102 i=1,N-1
       s_rhs_1(i,j)=asy(i,j)*s(i,j-1)
     >             +(bsy(i,j)-1.d0)*s(i,j)
     >             +csy(i,j)*s(i,j+1)
 102  continue
       s_rhs_1(N,j)=0.d0
 101  continue
        
      do 103 i=1,N-1
      do 103 j=1,N-1
       s_rhs_2(i,j)=asx(i,j)*s_rhs_1(i-1,j)
     >             +(bsx(i,j)-1.d0)*s_rhs_1(i,j)
     >             +csx(i,j)*s_rhs_1(i+1,j)
 103  continue
c     Added terms are calculated [ s_rhs_2(i,j) ]

c     Calculate total RHS [ ds(i,j) ]
      do 104 i=1,N-1
      do 104 j=1,N-1
       ds(i,j)=s(i,j)
     >        +dt*v(i,j)
     >        +s_rhs_2(i,j)
 104  continue

c     Solve for introduced variable ff(i,j)
      do 105 j=1,N-1

c     Note: ff(0,j)=0 and ff(N,j)=0

c     Forward elimination
      do 106 i=2,N-1
       bsx(i,j)=bsx(i,j)-asx(i,j)*csx(i-1,j)/bsx(i-1,j)
       ds(i,j)=ds(i,j)-asx(i,j)*ds(i-1,j)/bsx(i-1,j)
 106  continue

c     Substitute for the last point
       ff(N-1,j)=ds(N-1,j)/bsx(N-1,j)

c     Backward substitution
      do 107 i=N-2,1,-1
       ff(i,j)=(ds(i,j)-ff(i+1,j)*csx(i,j))/bsx(i,j)
 107  continue

 105  continue
c     Note: Now ff(i,j) becomes the RHS

c     Solve for streamfunction s(i,j)
      do 108 i=1,N-1

c     Note: s(i,0)=0 and s(i,N)=0

c     Forward elimination
      do 109 j=2,N-1
       bsy(i,j)=bsy(i,j)-asy(i,j)*csy(i,j-1)/bsy(i,j-1)
       ff(i,j)=ff(i,j)-asy(i,j)*ff(i,j-1)/bsy(i,j-1)
 109  continue

c     Substitute for the last point
       s(i,N-1)=ff(i,N-1)/bsy(i,N-1)

c     Backward substitution
      do 110 j=N-2,1,-1
       s(i,j)=(ff(i,j)-s(i,j+1)*csy(i,j))/bsy(i,j)
 110  continue

 108  continue

C     SOLVE THE VORTICITY EQUATION

C     Calculate vorticity at the wall
C     NOTE:For these boundary conditions please refer to:
C          T. Stortkuhl, C. Zenger, S. ZN-1mer, "An Asymptotic Solution for
C          the Singularity at the Angular Point of the Lid Driven Cavity",
C          International Journal of Numerical Methods for Heat & Fluid Flow
C          1994, Vol 4, pp 47--59
      do 30 k=1,N-1
       v(k,0)=(
     >       -(s(k-1,1)+s(k,1)+s(k+1,1))/(3.d0*dh**2.)
     >       -(0.5d0*v(k-1,0)+0.5d0*v(k+1,0)
     >       +0.25d0*v(k-1,1)+v(k,1)+0.25d0*v(k+1,1))/(9.d0)
     >        )*9.d0/2.d0
       v(k,N)=(-1.d0/dh
     >       -(s(k-1,N-1)+s(k,N-1)+s(k+1,N-1))/(3.d0*dh**2.)
     >       -(0.5d0*v(k-1,N)+0.5d0*v(k+1,N)
     >       +0.25d0*v(k-1,N-1)+v(k,N-1)+0.25d0*v(k+1,N-1))/(9.d0)
     >        )*9.d0/2.d0
       v(0,k)=(
     >       -(s(1,k-1)+s(1,k)+s(1,k+1))/(3.d0*dh**2.)
     >       -(0.5d0*v(0,k-1)+0.5d0*v(0,k+1)
     >       +0.25d0*v(1,k-1)+v(1,k)+0.25d0*v(1,k+1))/(9.d0)
     >        )*9.d0/2.d0
       v(N,k)=(
     >       -(s(N-1,k-1)+s(N-1,k)+s(N-1,k+1))/(3.d0*dh**2.)
     >       -(0.5d0*v(N,k-1)+0.5d0*v(N,k+1)
     >       +0.25d0*v(N-1,k-1)+v(N-1,k)+0.25d0*v(N-1,k+1))/(9.d0)
     >        )*9.d0/2.d0
 30   continue
       v(0,0)=(-(s(1,1))/(3.d0*dh**2.)-(0.5d0*v(1,0)+0.5d0*v(0,1)
     >       +0.25d0*v(1,1))/(9.d0))*9.d0
       v(N,0)=(-(s(N-1,1))/(3.d0*dh**2.)-(0.5d0*v(N-1,0)+0.5d0*v(N,1)
     >       +0.25d0*v(N-1,1))/(9.d0))*9.d0
       v(N,N)=(-0.5d0/dh-(s(N-1,N-1))/(3.d0*dh**2.)-(0.5d0*v(N-1,N)
     >       +0.5d0*v(N,N-1)+0.25d0*v(N-1,N-1))/(9.d0))*9.d0
       v(0,N)=(-0.5d0/dh-(s(1,N-1))/(3.d0*dh**2.)-(0.5d0*v(1,N)
     >       +0.5d0*v(0,N-1)+0.25d0*v(1,N-1))/(9.d0))*9.d0

c     Construction of matrices for vorticity equation
      do 300 j=1,N-1
       avy(0,j)=-dt/dh**2.
       bvy(0,j)=1.d0+dt*2.d0/dh**2.
       cvy(0,j)=-dt/dh**2.
      do 301 i=1,N-1
       sx=(s(i+1,j)-s(i-1,j))/(2.d0*dh)
       sy=(s(i,j+1)-s(i,j-1))/(2.d0*dh)
        
       avx(i,j)=-dt/dh**2.-dt*Re*sy/(2.d0*dh)
       bvx(i,j)=1.d0+dt*2.d0/dh**2.
       cvx(i,j)=-dt/dh**2.+dt*Re*sy/(2.d0*dh)

       avy(i,j)=-dt/dh**2.+dt*Re*sx/(2.d0*dh)
       bvy(i,j)=1.d0+dt*2.d0/dh**2.
       cvy(i,j)=-dt/dh**2.-dt*Re*sx/(2.d0*dh)
 301  continue
       avy(N,j)=-dt/dh**2.
       bvy(N,j)=1.d0+dt*2.d0/dh**2.
       cvy(N,j)=-dt/dh**2.
 300  continue

c     Calculate the added terms to RHS
      do 302 j=1,N-1
      do 302 i=0,N
       v_rhs_1(i,j)=avy(i,j)*v(i,j-1)
     >             +(bvy(i,j)-1.d0)*v(i,j)
     >             +cvy(i,j)*v(i,j+1)
 302  continue
        
      do 303 i=1,N-1
      do 303 j=1,N-1
       v_rhs_2(i,j)=avx(i,j)*v_rhs_1(i-1,j)
     >             +(bvx(i,j)-1.d0)*v_rhs_1(i,j)
     >             +cvx(i,j)*v_rhs_1(i+1,j)
 303  continue
c     Added terms are calculated [ v_rhs_2(i,j) ]

c     Calculate total RHS [ dv(i,j) ]
      do 304 i=1,N-1
      do 304 j=1,N-1
       dv(i,j)=v(i,j)
     >        +v_rhs_2(i,j)
 304  continue

        
c     Solve for introduced variable gg(i,j)
      do 305 j=1,N-1

c     Apply BC
       gg(0,j)=avy(0,j)*v(0,j-1)
     >        +bvy(0,j)*v(0,j)
     >        +cvy(0,j)*v(0,j+1)
       gg(N,j)=avy(N,j)*v(N,j-1)
     >        +bvy(N,j)*v(N,j)
     >        +cvy(N,j)*v(N,j+1)
       dv(1,j)=dv(1,j)-avx(1,j)*gg(0,j)
       dv(N-1,j)=dv(N-1,j)-cvx(N-1,j)*gg(N,j)

c     Forward elimination
      do 306 i=2,N-1
       bvx(i,j)=bvx(i,j)-avx(i,j)*cvx(i-1,j)/bvx(i-1,j)
       dv(i,j)=dv(i,j)-avx(i,j)*dv(i-1,j)/bvx(i-1,j)
 306  continue

c     Substitute for the last point
       gg(N-1,j)=dv(N-1,j)/bvx(N-1,j)

c     Backward substitution
      do 307 i=N-2,1,-1
       gg(i,j)=(dv(i,j)-gg(i+1,j)*cvx(i,j))/bvx(i,j)
 307  continue

 305  continue
c     Note: Now gg(i,j) becomes the RHS

c     Solve for vorticity v(i,j)
      do 308 i=1,N-1

c     Apply BC
       gg(i,1)=gg(i,1)-avy(i,1)*v(i,0)
       gg(i,N-1)=gg(i,N-1)-cvy(i,N-1)*v(i,N)

c     Forward elimination
      do 309 j=2,N-1
       bvy(i,j)=bvy(i,j)-avy(i,j)*cvy(i,j-1)/bvy(i,j-1)
       gg(i,j)=gg(i,j)-avy(i,j)*gg(i,j-1)/bvy(i,j-1)
 309  continue

c     Substitute for the last point
       v(i,N-1)=gg(i,N-1)/bvy(i,N-1)

c     Backward substitution
      do 310 j=N-2,1,-1
       v(i,j)=(gg(i,j)-v(i,j+1)*cvy(i,j))/bvy(i,j)
 310  continue

 308  continue

C     Check the residuals to see if convergence is achieved. 
C     You can comment out all or some, for a faster run.

c     residual_1 is the residual of the governing equations.
       residual_1_s_A=0.d0
       residual_1_v_A=0.d0
c     residual_2 is the change in variables (indicates the significant digit) in a time step.
       residual_2_s_A=0.d0
       residual_2_v_A=0.d0
c     residual_3 is the normalized change in variables (indicates percent change) in a time step.
       residual_3_s_A=0.d0
       residual_3_v_A=0.d0
      
      do 60 i=1,N-1
      do 60 j=1,N-1
       residual_1_s_B=abs(
     >   (s(i-1,j)-2.d0*s(i,j)+s(i+1,j))/dh**2.
     >  +(s(i,j-1)-2.d0*s(i,j)+s(i,j+1))/dh**2.
     >  +v(i,j)          )
       residual_1_v_B=abs(
     >   (1.d0/Re)*(v(i-1,j)-2.d0*v(i,j)+v(i+1,j))/dh**2.
     >  +(1.d0/Re)*(v(i,j-1)-2.d0*v(i,j)+v(i,j+1))/dh**2.
     >  -(s(i,j+1)-s(i,j-1))/(2.d0*dh)*(v(i+1,j)-v(i-1,j))/(2.d0*dh)
     >  +(s(i+1,j)-s(i-1,j))/(2.d0*dh)*(v(i,j+1)-v(i,j-1))/(2.d0*dh)
     >                   )

       residual_2_s_B=abs(s(i,j)-s_old(i,j))
       residual_2_v_B=abs(v(i,j)-v_old(i,j))

       residual_3_s_B=abs((s(i,j)-s_old(i,j))/s_old(i,j))
       residual_3_v_B=abs((v(i,j)-v_old(i,j))/v_old(i,j))

       residual_1_s_A=max(residual_1_s_A,residual_1_s_B)
       residual_1_v_A=max(residual_1_v_A,residual_1_v_B)

       residual_2_s_A=max(residual_2_s_A,residual_2_s_B)
       residual_2_v_A=max(residual_2_v_A,residual_2_v_B)

       residual_3_s_A=max(residual_3_s_A,residual_3_s_B)
       residual_3_v_A=max(residual_3_v_A,residual_3_v_B)
 60   continue
      
C     Output the residuals
       write(*,*) ' '
       write(*,*) iteration
       write(*,*) residual_1_s_A,residual_1_v_A
       write(*,*) residual_2_s_A,residual_2_v_A
       write(*,*) residual_3_s_A,residual_3_v_A
c     condition to stop iterations
       if((residual_1_s_A.lt.1.d-6).AND.(residual_1_v_A.lt.1.d-6))
     > goto 1000
      
 999  continue
      
C     Record the CPU time at finish
 1000 call cpu_time(time_finish)

       write(*,*) ' '
       write(*,*) 'Convergence is achieved in',iteration,'   iterations'
       write(*,*) 'CPU time=',time_finish-time_start

c     Output to a file        
      open(1,file='out.txt')
      do 1111 i=0,N
      do 1111 j=0,N
       write(1,2222) x(i),y(j),s(i,j),v(i,j)
 1111 continue
 2222 format(f8.4,x,f8.4,x,es25.18,x,es25.18)

      stop
      end