triacavity.f

C     This program solves the 2-D steady incompressible driven triangular cavity flow.
C     The Navier-Stokes equations in streamfunction and vorticity formulation
C     are solved using Successive Over Relaxation SOR method.
C     Visit http://www.cavityflow.com
C
C     s(i,j) ==>  streamfunction variable
C     v(i,j) ==>  vorticity variable
C     x(i,j) ==>  x-coordinate
C     y(j,j) ==>  y-coordinate
C     dh     ==>  grid spacing
C     Re     ==>  Reynolds Number
C     ax,bx,h,cx,cy  ==>  corner points of triangle
      
      
      program main
      
      implicit double precision (a-h,o-z)
      
      parameter(N=512)
      
      common / flow variables /
     >s(0:N,0:N),v(0:N,0:N)
      
      common / geometry /
     >xx(0:N),yy(0:N),x(0:N,0:N),y(0:N,0:N)
      
      do 1 i=0,N
      do 1 j=0,N
      s(i,j)=0.d0
      v(i,j)=0.d0
 1    continue
      
      
      Re=500.d0
      
      beta=1.0d0
      
      
C     Read data from this file as an initial guess
      open(1,file='triangle_1_0200.bin',form="unformatted",
     >access="stream")
      read(1) ax
      read(1) bx
      read(1) h
      read(1) cx
      read(1) cy
      read(1) Redummy
      do 2 i=0,N
      do 2 j=0,N-i
      read(1) s(i,j),v(i,j)
 2    continue
      close(1)
      
      dzdx=1.d0/(bx-ax)
      dzdy=(ax-cx)/((bx-ax)*(cy-h))
      dndy=1.d0/(cy-h)
      
      dh=1.d0/dble(N)
      
      do 3 k=0,N
      xx(k)=dble(k)/dble(N)
      yy(k)=dble(k)/dble(N)
 3    continue
      
      do 4 i=0,N
      do 4 j=0,N
      y(i,j)=(yy(j)+h/(cy-h))/dndy
      x(i,j)=(xx(i)-dzdy*y(i,j)-(h*cx-ax*cy)/((bx-ax)*(cy-h)))/dzdx
 4    continue
      
      A=dzdx**2.+dzdy**2.
      B=dndy**2.
      C=2.d0*dzdy*dndy
      D=dzdx*dndy
      
      
      do 999 iteration=1,1000
      do 99 iter=1,25000
      
C     Iterate on Streamfunction
      do 5 j=1,N-2
      do 5 i=1,N-1-j
      s(i,j)=beta*(
     >( A*(s(i-1,j)+s(i+1,j))
     > +B*(s(i,j-1)+s(i,j+1))
     > +0.25d0*C*(s(i+1,j+1)+s(i-1,j-1)-s(i-1,j+1)-s(i+1,j-1))
     > +dh**2.*v(i,j)
     >)/(2.d0*A+2.d0*B)
     >            )+(1.d0-beta)*s(i,j)
 5    continue
      
C     Calculate Vorticity at Boundaries
      do 6 i=1,N-1
      v(i,0)=-B*(2.d0*s(i,1)-2.d0*dh/dndy)/dh**2.
      v(0,i)=-A*2.d0*s(1,i)/dh**2.
 6    continue
      do 7 i=1,N-1
      v(i,N-i)=-1.d0*(2.d0*A*s(i-1,N-i)
     >               +2.d0*B*s(i,N-i-1)
     >               +2.d0*0.25d0*C*s(i-1,N-i-1))/dh**2.
 7    continue
      
C     Apply perfect slip condition at the corner points
      v(0,0)=0.0d0
      v(N,0)=0.0d0
      v(0,N)=0.0d0
      
C     Iterate on Vorticity
      do 8 j=1,N-2
      do 8 i=1,N-1-j
      v(i,j)=beta*(
     >( A*(v(i-1,j)+v(i+1,j))
     > +B*(v(i,j-1)+v(i,j+1))
     > +0.25d0*C*(v(i+1,j+1)+v(i-1,j-1)-v(i-1,j+1)-v(i+1,j-1))
     > -Re*D*((s(i,j+1)-s(i,j-1))/2.d0)*((v(i+1,j)-v(i-1,j))/2.d0)
     > +Re*D*((s(i+1,j)-s(i-1,j))/2.d0)*((v(i,j+1)-v(i,j-1))/2.d0)
     >)/(2.d0*A+2.d0*B)
     >            )+(1.d0-beta)*v(i,j)
 8    continue
      
 99   continue
      
      res_s_A=0.d0
      res_v_A=0.d0
      
      do 9 j=1,N-2
      do 9 i=1,N-1-j
      
      res_s_B=abs(
     > A*(s(i-1,j)-2.d0*s(i,j)+s(i+1,j))/dh**2.
     >+B*(s(i,j-1)-2.d0*s(i,j)+s(i,j+1))/dh**2.
     >+0.25d0*C*(s(i+1,j+1)+s(i-1,j-1)-s(i-1,j+1)-s(i+1,j-1))/dh**2.
     >+v(i,j)
     >           )
      res_v_B=abs(
     > A*(v(i-1,j)-2.d0*v(i,j)+v(i+1,j))/(Re*dh**2.)
     >+B*(v(i,j-1)-2.d0*v(i,j)+v(i,j+1))/(Re*dh**2.)
     >+0.25d0*C*(v(i+1,j+1)+v(i-1,j-1)
     >          -v(i-1,j+1)-v(i+1,j-1))/(Re*dh**2.)
     >-D*((s(i,j+1)-s(i,j-1))/(2.0d0*dh))
     >           *((v(i+1,j)-v(i-1,j))/(2.0d0*dh))
     >+D*((s(i+1,j)-s(i-1,j))/(2.0d0*dh))
     >           *((v(i,j+1)-v(i,j-1))/(2.0d0*dh))
     >           )
      
      res_s_A=max(res_s_A,res_s_B)
      res_v_A=max(res_v_A,res_v_B)
      
 9    continue
      
C     Output the residuals
      open(2,file='residuals.txt',ACCESS='APPEND')
      write(2,*) 'iteration number=',iteration*50000
      write(2,*) 'res_f=',res_s_A,'   res_g=',res_v_A
      close(2)
      
C     Output the variables
      open(3,file='out.dat')
      rewind(3)
	do 10 i=0,N
      do 10 j=0,N
      write(3,80) x(i,j),y(i,j),s(i,j),v(i,j)
 10   continue
      close(3)
 80   format(f8.4,x,f8.4,x,es25.18,x,es25.18)
      
C     Condition to stop the code
      if((res1s_A.lt.1.d-9).AND.(res1v_A.lt.1.d-9)) then
      goto 100
      endif
      
 999  continue
      
 100  continue
      stop
      end