Wednesday, December 3, 2008

Velocity potential by Relaxation Tecnique in fortran

c Program to find Velocity potential by Relaxation Tecnique

dimension u(10,10),v(10,10),div(10,10),ngrid(10,10),ogrid(10,10)
real forcef,ngrid,ogrid,residue
integer n,k

print*,'Enter the size of grid excluding boundary:'
read*,n
print*,'Enter the distance between grid points :'
read*,d


do i=1,n+2
u(1,i)=0
u(i,1)=0
u(n+2,i)=0
u(i,n+2)=0
v(1,i)=0
v(i,1)=0
v(n+2,i)=0
v(i,n+2)=0
div(1,i)=0
div(i,1)=0
div(n+2,i)=0
div(i,n+2)=0
ngrid(1,i)=0
ngrid(i,1)=0
ngrid(n+2,i)=0
ngrid(i,n+2)=0
ogrid(1,i)=0
ogrid(i,1)=0
ogrid(n+2,i)=0
ogrid(i,n+2)=0
enddo

open (1,file='uwnd.dat')
read (1,*)((u(i,j),j=2,n+1),i=2,n+1)
close(1)
open (2,file='vwnd.dat')
read (2,*)((v(i,j),j=2,n+1),i=2,n+1)
close(2)
do i=2,n+1
do j=2,n+1
div(i,j)=((u(i+1,j)-u(i-1,j))+(v(i,j+1)-v(i,j-1)))/(2*d)
enddo
enddo

write(*,*)'Divergence values'
do i=2,n+1
write(*,*) (div(i,j),j=2,n+1)
enddo

do i=2,n+1
do j=2,n+1
ngrid(i,j)=0
enddo
enddo

k=1
10 print*,'Trial :',k

do i=2,n+1
do j=2,n+1
ogrid(i,j)=ngrid(i,j)
enddo
enddo

do i=2,n+1
do j=2,n+1
residue=ogrid(i-1,j)+ogrid(i+1,j)+ogrid(i,j-1)+ogrid(i,j+1)-(4*ogrid(i,j))+(d*d*div(i,j))

ngrid(i,j)=ogrid(i,j)+(residue/4)
enddo
enddo

do i=2,n+1
do j=2,n+1
if((ngrid(i,j)-ogrid(i,j)).gt.0.0001) then
k=k+1
goto 10
endif
enddo
enddo

open(3,file="velocity_potential.txt")
do i=2,n+1
write(3,*)(ngrid(i,j),j=2,n+1)
enddo
write(*,*)'Velocity potentials are written in the file velocity_potential.txt'
do i=2,n+1
write(*,*) (ngrid(i,j),j=2,n+1)
enddo
stop
end

No comments: