sorcavity.f

C****************************************************************************C
C                                                                            C
C     This program solve the Navier-Stokes (NS) equations using the          C
C     Successive Over Relaxation (SOR) method. The solution is spatially     C
C     second order accurate. 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     ==> Reynolds Number             C
C                                            C
C********************************************C

      
      program main
      
      implicit double precision (a-h,o-z)
      
      parameter(N=1024)
      
      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)
      
      
      Re=2500.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     Use the beta values given in Table 1
      beta=1.0d0
      
C     Read data from this file as an initial guess
      open(1,file='sor_Re01000.bin',form="unformatted",access="stream")
C     Read a dummy number (each file contains the Reynolds number as the first data)
      read(1) redummy
      do 2 i=0,N
      do 2 j=0,N
      read(1) s(i,j),v(i,j)
 2    continue
      close(1)
      
      
      do 999 iteration=1,1000
      do 99 iter=1,50000
      
      do 3 i=0,N
      do 3 j=0,N
      s_old(i,j)=s(i,j)
      v_old(i,j)=v(i,j)
 3    continue
      
C     Iterate on Streamfunction
      do 4 i=1,N-1
      do 4 j=1,N-1
      s(i,j)=beta*(0.25d0*(
     > s(i-1,j)+s(i+1,j)+s(i,j-1)+s(i,j+1)
     >+dh**2.*v(i,j)
     >                    )
     >            )+(1.d0-beta)*s(i,j)
 4    continue
      
C     Calculate Vorticity at Boundaries using Jensen's formula
      do 5 i=1,N-1
      v(i,0)=0.5d0*(-8.d0*s(i,1)+s(i,2))/dh**2.
      v(i,N)=0.5d0*(-8.d0*s(i,N-1)+s(i,N-2))/dh**2.-3.d0/dh
 5    continue
      do 6 j=1,N-1
      v(0,j)=0.5d0*(-8.d0*s(1,j)+s(2,j))/dh**2.
      v(N,j)=0.5d0*(-8.d0*s(N-1,j)+s(N-2,j))/dh**2.
 6    continue
      
C     Calculate Vorticity at the corners
C     NOTE:We do not need these points for interior solution
C          We only need these points for plotting purposes
C          For these boundary conditions please refer to:
C          T. Stortkuhl, C. Zenger, S. Zimmer, "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
      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     Iterate on Vorticity
      do 7 i=1,N-1
      do 7 j=1,N-1
      v(i,j)=beta*(0.25d0*(
     > v(i-1,j)+v(i+1,j)+v(i,j-1)+v(i,j+1)
     >-0.25d0*Re*(s(i,j+1)-s(i,j-1))*(v(i+1,j)-v(i-1,j))
     >+0.25d0*Re*(s(i+1,j)-s(i-1,j))*(v(i,j+1)-v(i,j-1))
     >                    )
     >            )+(1.d0-beta)*v(i,j)
 7    continue
      
 99   continue
      
      residual_1_s_A=0.d0
      residual_1_v_A=0.d0
      
      residual_2_s_A=0.d0
      residual_2_v_A=0.d0
      
      residual_3_s_A=0.d0
      residual_3_v_A=0.d0
           
      do 8 i=1,N-1
      do 8 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_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_B=abs(s(i,j)-s_old(i,j))
      residual_2_v_B=abs(v(i,j)-v_old(i,j))
      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_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_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)
 8    continue
      
C     Output the residuals
      open(2,file='residuals.txt',ACCESS='APPEND')
      write(2,*) 'iteration number=',iteration*50000
      write(2,*) 'residual_1_s=',residual_1_s_A,
     >        '   residual_1_v=',residual_1_v_A
      write(2,*) 'residual_2_s=',residual_2_s_A,
     >        '   residual_2_v=',residual_2_v_A
      write(2,*) 'residual_3_s=',residual_3_s_A,
     >        '   residual_3_v=',residual_3_v_A
      close(2)
      
C     Output the variables
      open(3,file='out.txt')
      rewind(3)
      do 9 i=0,N
      do 9 j=0,N
      write(3,80) x(i),y(j),s(i,j),v(i,j)
 9    continue
      close(3)
 80   format(f8.4,x,f8.4,x,es25.18,x,es25.18)
      
C     Condition to stop the code
      if((residual_1_s_A.lt.1.d-9).AND.(residual_1_v_A.lt.1.d-9))goto 10
      
      
 999  continue
      
 10   continue
      stop
      end