Differential mesh (FORTRAN VERSION) -- CFD

[ Show ]

Support VoyForums

[ Shrink ]

VoyForums Announcement: Programming and providing support for this service has been a labor of love since 1997. We are one of the few services online who values our users' privacy, and have never sold your information. We have even fought hard to defend your privacy in legal cases; however, we've done it with almost no financial support -- paying out of pocket to continue providing the service. Due to the issues imposed on us by advertisers, we also stopped hosting most ads on the forums many years ago. We hope you appreciate our efforts.

Show your support by donating any amount. (Note: We are still technically a for-profit company, so your contribution is not tax-deductible.) PayPal Acct: Feedback:

Donate to VoyForums (PayPal):

[ Login ] [ Contact Forum Admin ] [ Main index ] [ Post a new message ] [ Search | Check update time ]

[ Next Thread | Previous Thread | Next Message | Previous Message ]

Date Posted: 08:53:39 05/18/04 Tue
Author: Reza Khaksar Fard
Subject: Differential mesh (FORTRAN VERSION)

write(*,*) 'differential mesh ? (type 1 for yes)'
read (*,*) ldif
if (ldif.eq.1)
& CALL MESH (xlap,ylap, nxcells-3, Nycells-1)

c unstretched .
SUBROUTINE MESH(x,y,Ni,Nj)

dimension x(100,100),y(100,100)
dimension xnew(100,100),ynew(100,100),gx(100),gy(100),gama(100)
dimension a(100),b(100),c(100),d(100),alpha(100),beta(100)
eps=.001
omega=1.0

do 91 j=1,Nj
do 92 i=1,Ni+2
xnew(i,j)=x(i,j)
ynew(i,j)=y(i,j)
92 continue
91 continue

c ****************************************
c The differential grid (Unstretched) :

iterno=0
70 DO 61 j=2,Nj-1
DO 11 i=2,Ni+1
xeta=D1(x(i,j+1),x(i,j-1),1.)
yeta=D1(y(i,j+1),y(i,j-1),1.)
xsi=D1(x(i+1,j),x(i-1,j),1.)
ysi=D1(y(i+1,j),y(i-1,j),1.)
alpha(i-1)=xeta**2+yeta**2
beta(i-1)=xsi*xeta+ysi*yeta
gama(i-1)=xsi**2+ysi**2
gx(i-1)=beta(i-1)*2.*D3(x(i+1,j+1),x(i+1,j-1),x(i-1,j+1),
& x(i-1,j-1),1.,1.)-gama(i-1)*(x(i,j+1)+x(i,j-1))
gy(i-1)=beta(i-1)*2.*D3(y(i+1,j+1),y(i+1,j-1),y(i-1,j+1),
& y(i-1,j-1),1.,1.)-gama(i-1)*(y(i,j+1)+y(i,j-1))
11 CONTINUE

DO 21 k=1,Ni
d(k)=-2.*(alpha(k)+gama(k))
a(k)=alpha(k)
b(k)=alpha(k)
c(k)=gx(k)
21 CONTINUE
c(1)=c(1)-alpha(1)*x(1,j)
c(Ni)=c(Ni)-alpha(Ni)*x(Ni+2,j)
CALL SY(1,Ni,b,d,a,c)
DO 31 k=1,Ni
xnew(k+1,j)=c(k)
31 CONTINUE

DO 41 k=1,Ni
d(k)=-2.*(alpha(k)+gama(k))
c(k)=gy(k)
41 CONTINUE
c(1)=c(1)-alpha(1)*y(1,j)
c(Ni)=c(Ni)-alpha(Ni)*y(Ni+2,j)
CALL SY(1,Ni,b,d,a,c)
DO 51 k=1,Ni
ynew(k+1,j)=c(k)
51 CONTINUE
61 CONTINUE

c The comparison :

DeltaX=0.
DO 81 i=1,Ni+2
DO 71 j=2,Nj-1
c deltax=ABS(x(i,j)-xnew(i,j))
c deltay=ABS(y(i,j)-ynew(i,j))
c IF (deltax.GT.eps) GO TO 12
c IF (deltay.GT.eps) GO TO 12
DeltaX=DeltaX+(ABS(x(i,j)-xnew(i,j))+ ABS(y(i,j)-ynew(i,j)))
c IF (deltax.GT.eps) GO TO 12
c IF (deltay.GT.eps) GO TO 12
71 CONTINUE
81 CONTINUE
DeltaX=log10(DeltaX/(Ni+2)/(Nj-2))
if(DeltaX .gt. -3.) GO TO 12

GO TO 100

c The S.L.O.R. :

c 12 write(*,*) 'iteration # ',iterno,deltax,deltay,i,j
12 write(*,*) 'iteration # ',iterno,deltax
DO 32 i=1,Ni+2
DO 22 j=2,Nj-1
x(i,j)=x(i,j)+omega*(xnew(i,j)-x(i,j))
y(i,j)=y(i,j)+omega*(ynew(i,j)-y(i,j))
22 CONTINUE
32 CONTINUE

iterno=iterno+1

GO TO 70

100 RETURN
END

c ########################################################################

SUBROUTINE SY(il,iu,bb,dd,aa,cc)
dimension aa(1),bb(1),cc(1),dd(1)
lp=il+1
DO 10 i=lp,iu
r=bb(i)/dd(i-1)
dd(i)=dd(i)-r*aa(i-1)
cc(i)=cc(i)-r*cc(i-1)
10 CONTINUE
cc(iu)=cc(iu)/dd(iu)
DO 20 i=lp,iu
j=iu-i+il
cc(j)=(cc(j)-aa(j)*cc(j+1))/dd(j)
20 CONTINUE
RETURN
END

FUNCTION D1(yplus,yminus,def)
c Central difference for a 1st order derivative term :
D1=(yplus-yminus)/2./def
RETURN
END

FUNCTION D3(ypp,ypm,ymp,ymm,def1,def2)
c Central difference for a 2st order mixed derivative term :
D3=(ypp-ypm-ymp+ymm)/4./def1/def2
RETURN
END

[ Next Thread | Previous Thread | Next Message | Previous Message ]