c---------------------------------------------------------------------- c This is a general purpose time marching 2-D Navier-Stokes solver c c input data must be in a file 'chan.dat' in the run directory c c program allows changes to occur during run time to c Reynolds number, the program can be instructed to c save everything and then be restarted later. c c Uses multigrid to solve for stream function c from vorticity and different methods can be used c for time marching the vorticity c c-------------------------------------------------------------channel-- program channel include 'step.h' call start 1 continue call runtimesequence 2 call whatnext(key) if(key.gt.0)then if(key.eq.5)then goto 2 else goto 1 endif endif call shutdown stop end c-------------------------------------------------------------whatnext- subroutine whatnext(key) include 'step.h' write(logfile,22) 22 format(' 0=end,1=reflect,2=change Re,3=set dt,4=method,5=archive') read(inputfile,*)key if (key.eq.swapflow)then write(logfile,23) 23 format('swapping:') call swap else if(key.eq.changereynolds)then read(inputfile,*)reynolds write(logfile,24)reynolds 24 format('reynolds changed to:',f10.2) pex=dt/(reynolds*strouhal*dx*dx) pey=dt/(reynolds*strouhal*dy*dy) else if(key.eq.changedt)then read(inputfile,*)newdt write(logfile,25)newdt 25 format('dt changed to:',f10.4) dtratio=dt/newdt printsteps=dtratio*printsteps dtratio=newdt/dt do 30 i=0,ix1 do 30 j=-iy1,channeltype*iy1 oldvort(i,j)=dtratio*oldvort(i,j)+(1.d0-dtratio)*vort(i,j) 30 continue dt=newdt cx=dt/(strouhal*dx) cy=dt/(strouhal*dy) pex=dt/(reynolds*strouhal*dx*dx) pey=dt/(reynolds*strouhal*dy*dy) else if(key.eq.changemethod)then read(inputfile,*)method write(logfile,35)method 35 format('method changed to:',i4) else if(key.eq.archiveflow)then call archive write(logfile,45)reynolds 45 format('writing to archive re=',f10.2) endif return end c---------------------------------------------------------------------- subroutine swap include 'step.h' c c reflect y values c do 10 i=0,ix1 do 10 j=-iy1,iy1/2-1 temp=psi(i,j) psi(i,j)=-psi(i,iy1-j) psi(i,iy1-j)=-temp temp=newvort(i,j) newvort(i,j)=-newvort(i,iy1-j) newvort(i,iy1-j)=-temp temp=vort(i,j) vort(i,j)=-vort(i,iy1-j) vort(i,iy1-j)=-temp temp=oldvort(i,j) oldvort(i,j)=-oldvort(i,iy1-j) oldvort(i,iy1-j)=-temp temp=u(i,j) u(i,j)=u(i,iy1-j) u(i,iy1-j)=temp temp=v(i,j) v(i,j)=-v(i,iy1-j) v(i,iy1-j)=-temp 10 continue return end c----------------------------------------------------------------------- subroutine shutdown include 'step.h' write(psifile,30)ix1/2+1,(1+channeltype)*iy1/2+1 write(vortfile,30)ix1/2+1,iy1+1 30 format('aspect:.2',/,'height:0. -1. 1.',/,'contour ',2i4) write(psifile,31)((psi(i,j),i=0,ix1,2),j=-iy1,channeltype*iy1,2) write(vortfile,31)((vort(i,j),i=0,ix1,2),j=-iy1,channeltype*iy1,2) 31 format(8e14.5) close(unit=psifile) close(unit=vortfile) close(unit=sumfile) close(unit=logfile) close(unit=inputfile) return end c----------------------------------------------------------------------- subroutine runtimesequence include 'step.h' do 10 k=1,ns1 call onetimestep call export(1) if((k.gt.10).and.(vincmax.lt.1.e-6))then write(logfile,5)k,reynolds 5 format('converged to steady flow after ',i5,'steps : Re =',f8.2) goto 15 endif 10 continue 15 do 20 k=1,ns2 call onetimestep call export(2) 20 continue return end c----------------------------------------------------------------------- subroutine onetimestep include 'step.h' currenttime=currenttime+dt call setflux if (method.eq.central)then call stepcentral else if(method.eq.upwind)then call stepupwind else if(method.eq.ncupwind)then call stepnonc else if(method.eq.laxwendroff)then call steplax else if(method.eq.quickest)then call stepquick endif call calculatepsi call calculatevelocities call calculateboundaryvorticity do 50 i=0,ix1 do 50 j=-iy1,channeltype*iy1 oldvort(i,j)=vort(i,j) vort(i,j)=newvort(i,j) 50 continue return end c-------------------------------------------------------------------- subroutine export(sequence) include 'step.h' integer sequence s1=0. s2=0. do 10 i=ixc1,ix1-1 if(vort(i,-iy1)*vort(i+1,-iy1).lt.0.)then fac=(0.-vort(i,-iy1))/(vort(i+1,-iy1)-vort(i,-iy1)) s1=(i+fac-ixc1)*dx endif if(vort(i,channeltype*iy1)*vort(i+1,channeltype*iy1).lt.0.)then fac=(0.-vort(i,channeltype*iy1))/ + (vort(i+1,channeltype*iy1)-vort(i,channeltype*iy1)) s2=(i+fac-ixc1)*dx endif 10 continue pmax=0. pmin=0. ipmax=0 ipmin=0 jpmin=0 jpmax=0 do 12 i=0,ix1 do 12 j=-iy1,channeltype*iy1 if(psi(i,j).gt.pmax)then pmax=psi(i,j) ipmax=i jpmax=j endif if(psi(i,j).lt.pmin)then pmin=psi(i,j) ipmin=i jpmin=j endif 12 continue write(screen,15)currenttime,psi(0,iy1),reynolds, + vort(ixc1,0),s1,ipmin,jpmin,pmin, + ipmax,jpmax,pmax,vincmax 15 format(f8.4,f8.3,f8.2,f8.1,f8.3,2(2i4,f10.5),e14.4) if(sequence.eq.2)then call calculateasymmetry(asym,vasym) write(logfile,17)currenttime,reynolds, + s1,ipmin,jpmin,pmin, + s2,ipmax,jpmax,pmax write(sumfile,16)reynolds,s1,s2,asym,vasym 16 format(f8.2,2f10.4,2e14.4) 17 format(f8.3,f8.2,2(f8.4,2i4,f10.5)) endif return end c------------------------------------------------------------------ subroutine calculateasymmetry(asym,vasym) include 'step.h' area=ixc1*dx+3*(ix1-ixc1)*dx asym=0. do 20 i=1,ix1 do 10 j=iy1/2,channeltype*iy1 f=psi(i,j)+psi(i,iy1-j) asym=asym+f*f 10 continue 20 continue asym=sqrt(asym*dx*dy/area) vasym=(1.-vfac)*v(jv,iy1/2)+vfac*v(jv+1,iy1/2) asym=sign(asym,vasym) return end c--------------------------------------------------------------------- subroutine initialiseall() include 'step.h' do 20 i=0,ix1 do 10 j=-iy1,channeltype*iy1 psi(i,j)=0. vort(i,j)=0. oldvort(i,j)=0. newvort(i,j)=0. u(i,j)=0. v(i,j)=0. 10 continue 20 continue return end c---------------------------------------------------------------------- subroutine steplax include 'step.h' real *8 peu,pev vincmax=0. do 10 i=1,ixc1 do 10 j=1,iy1-1 au=-0.5*cx*u(i-1,j) bu=0. cu=0.5*cx*u(i+1,j) av=-0.5*cy*v(i,j-1) bv=0. cv=0.5*cy*v(i,j+1) peu=pex + cx*cx*u(i,j)*u(i,j)/2 pev=pey + cy*cy*v(i,j)*v(i,j)/2 vincrement(i,j)= + -au*vort(i-1,j)-cu*vort(i+1,j) + -av*vort(i,j-1)-cv*vort(i,j+1) + +peu*(vort(i+1,j)-2*vort(i,j)+vort(i-1,j)) + +pev*(vort(i,j-1)-2*vort(i,j)+vort(i,j+1)) newvort(i,j)=vort(i,j)+vincrement(i,j) temp=abs(vincrement(i,j)) if(temp.gt.vincmax)vincmax=temp 10 continue do 20 i=ixc1+1,ix1-1 do 20 j=-iy1+1,channeltype*iy1-1 au=-0.5*cx*u(i-1,j) bu=0. cu=0.5*cx*u(i+1,j) av=-0.5*cy*v(i,j-1) bv=0. cv=0.5*cy*v(i,j+1) peu=pex + cx*cx*u(i,j)*u(i,j)/2 pev=pey + cy*cy*v(i,j)*v(i,j)/2 vincrement(i,j)= + -au*vort(i-1,j)-cu*vort(i+1,j) + -av*vort(i,j-1)-cv*vort(i,j+1) + +peu*(vort(i+1,j)-2*vort(i,j)+vort(i-1,j)) + +pev*(vort(i,j-1)-2*vort(i,j)+vort(i,j+1)) newvort(i,j)=vort(i,j)+vincrement(i,j) temp=abs(vincrement(i,j)) if(temp.gt.vincmax)vincmax=temp 20 continue i=ix1 do 30 j=-iy1+1,channeltype*iy1-1 au=-0.5*cx*u(i-1,j) bu=0. cu=0.5*cx*u(i-1,j) av=-0.5*cy*v(i,j-1) bv=0. cv=0.5*cy*v(i,j+1) peu=pex + cx*cx*u(i,j)*u(i,j)/2 pev=pey + cy*cy*v(i,j)*v(i,j)/2 vincrement(i,j)= + -au*vort(i-1,j)-cu*vort(i-1,j) + -av*vort(i,j-1)-cv*vort(i,j+1) + +peu*(vort(i-1,j)-2*vort(i,j)+vort(i-1,j)) + +pev*(vort(i,j-1)-2*vort(i,j)+vort(i,j+1)) newvort(i,j)=vort(i,j)+vincrement(i,j) temp=abs(vincrement(i,j)) if(temp.gt.vincmax)vincmax=temp 30 continue return end c---------------------------------------------------------------------- subroutine stepcentral include 'step.h' vincmax=0. do 10 i=1,ixc1 do 10 j=1,iy1-1 au=-0.5*cx*u(i-1,j) bu=0. cu=0.5*cx*u(i+1,j) av=-0.5*cy*v(i,j-1) bv=0. cv=0.5*cy*v(i,j+1) vincrement(i,j)= + -au*vort(i-1,j)-cu*vort(i+1,j) + -av*vort(i,j-1)-cv*vort(i,j+1) + +pex*(vort(i+1,j)-2*vort(i,j)+vort(i-1,j)) + +pey*(vort(i,j-1)-2*vort(i,j)+vort(i,j+1)) newvort(i,j)=vort(i,j)+vincrement(i,j) temp=abs(vincrement(i,j)) if(temp.gt.vincmax)vincmax=temp 10 continue do 20 i=ixc1+1,ix1-1 do 20 j=-iy1+1,channeltype*iy1-1 au=-0.5*cx*u(i-1,j) bu=0. cu=0.5*cx*u(i+1,j) av=-0.5*cy*v(i,j-1) bv=0. cv=0.5*cy*v(i,j+1) vincrement(i,j)= + -au*vort(i-1,j)-cu*vort(i+1,j) + -av*vort(i,j-1)-cv*vort(i,j+1) + +pex*(vort(i+1,j)-2*vort(i,j)+vort(i-1,j)) + +pey*(vort(i,j-1)-2*vort(i,j)+vort(i,j+1)) newvort(i,j)=vort(i,j)+vincrement(i,j) temp=abs(vincrement(i,j)) if(temp.gt.vincmax)vincmax=temp 20 continue i=ix1 do 30 j=-iy1+1,channeltype*iy1-1 au=-0.5*cx*u(i-1,j) bu=0. cu=0.5*cx*u(i-1,j) av=-0.5*cy*v(i,j-1) bv=0. cv=0.5*cy*v(i,j+1) vincrement(i,j)= + -au*vort(i-1,j)-cu*vort(i-1,j) + -av*vort(i,j-1)-cv*vort(i,j+1) + +pex*(vort(i-1,j)-2*vort(i,j)+vort(i-1,j)) + +pey*(vort(i,j-1)-2*vort(i,j)+vort(i,j+1)) newvort(i,j)=vort(i,j)+vincrement(i,j) temp=abs(vincrement(i,j)) if(temp.gt.vincmax)vincmax=temp 30 continue return end c---------------------------------------------------------------------- subroutine stepnonc include 'step.h' vincmax=0. do 10 i=1,ixc1 do 10 j=1,iy1-1 if(u(i,j).ge.0.)then au=-cx*u(i,j) bu= cx*u(i,j) cu= 0. else au=0 . bu=-cx*u(i,j) cu= cx*u(i,j) endif if(v(i,j).ge.0.)then av=-cy*v(i,j) bv= cy*v(i,j) cv= 0. else av= 0. bv=-cy*v(i,j) cv= cy*v(i,j) endif vincrement(i,j)= + -au*vort(i-1,j)-(bu+bv)*vort(i,j)-cu*vort(i+1,j) + -av*vort(i,j-1)-cv*vort(i,j+1) + +pex*(vort(i+1,j)-2*vort(i,j)+vort(i-1,j)) + +pey*(vort(i,j-1)-2*vort(i,j)+vort(i,j+1)) newvort(i,j)=vort(i,j)+vincrement(i,j) temp=abs(vincrement(i,j)) if(temp.gt.vincmax)vincmax=temp 10 continue do 20 i=ixc1+1,ix1-1 do 20 j=-iy1+1,channeltype*iy1-1 if(u(i,j).ge.0.)then au=-cx*u(i,j) bu= cx*u(i,j) cu= 0. else au=0 . bu=-cx*u(i,j) cu= cx*u(i,j) endif if(v(i,j).ge.0.)then av=-cy*v(i,j) bv= cy*v(i,j) cv= 0. else av= 0. bv=-cy*v(i,j) cv= cy*v(i,j) endif vincrement(i,j)= + -au*vort(i-1,j)-(bu+bv)*vort(i,j)-cu*vort(i+1,j) + -av*vort(i,j-1)-cv*vort(i,j+1) + +pex*(vort(i+1,j)-2*vort(i,j)+vort(i-1,j)) + +pey*(vort(i,j-1)-2*vort(i,j)+vort(i,j+1)) newvort(i,j)=vort(i,j)+vincrement(i,j) temp=abs(vincrement(i,j)) if(temp.gt.vincmax)vincmax=temp 20 continue i=ix1 do 30 j=-iy1+1,channeltype*iy1-1 if(u(i,j).ge.0.)then au=-cx*u(i,j) bu= cx*u(i,j) cu= 0. else au=0 . bu=-cx*u(i,j) cu= cx*u(i,j) endif if(v(i,j).ge.0.)then av=-cy*v(i,j) bv= cy*v(i,j) cv= 0. else av= 0. bv=-cy*v(i,j) cv= cy*v(i,j) endif vincrement(i,j)= + -au*vort(i-1,j)-(bu+bv)*vort(i,j)-cu*vort(i-1,j) + -av*vort(i,j-1)-cv*vort(i,j+1) + +pex*(vort(i-1,j)-2*vort(i,j)+vort(i-1,j)) + +pey*(vort(i,j-1)-2*vort(i,j)+vort(i,j+1)) newvort(i,j)=vort(i,j)+vincrement(i,j) temp=abs(vincrement(i,j)) if(temp.gt.vincmax)vincmax=temp 30 continue return end c---------------------------------------------------------------------- subroutine stepupwind include 'step.h' vincmax=0. do 10 i=1,ixc1 do 10 j=1,iy1-1 if(u(i,j).ge.0.)then au=-cx*u(i-1,j) bu= cx*u(i,j) cu= 0. else au=0 . bu=-cx*u(i,j) cu= cx*u(i+1,j) endif if(v(i,j).ge.0.)then av=-cy*v(i,j-1) bv= cy*v(i,j) cv= 0. else av= 0. bv=-cy*v(i,j) cv= cy*v(i,j+1) endif vincrement(i,j)= + -au*vort(i-1,j)-(bu+bv)*vort(i,j)-cu*vort(i+1,j) + -av*vort(i,j-1)-cv*vort(i,j+1) + +pex*(vort(i+1,j)-2*vort(i,j)+vort(i-1,j)) + +pey*(vort(i,j-1)-2*vort(i,j)+vort(i,j+1)) newvort(i,j)=vort(i,j)+vincrement(i,j) temp=abs(vincrement(i,j)) if(temp.gt.vincmax)vincmax=temp 10 continue do 20 i=ixc1+1,ix1-1 do 20 j=-iy1+1,channeltype*iy1-1 if(u(i,j).ge.0.)then au=-cx*u(i-1,j) bu= cx*u(i,j) cu= 0. else au=0 . bu=-cx*u(i,j) cu= cx*u(i+1,j) endif if(v(i,j).ge.0.)then av=-cy*v(i,j-1) bv= cy*v(i,j) cv= 0. else av= 0. bv=-cy*v(i,j) cv= cy*v(i,j+1) endif vincrement(i,j)= + -au*vort(i-1,j)-(bu+bv)*vort(i,j)-cu*vort(i+1,j) + -av*vort(i,j-1)-cv*vort(i,j+1) + +pex*(vort(i+1,j)-2*vort(i,j)+vort(i-1,j)) + +pey*(vort(i,j-1)-2*vort(i,j)+vort(i,j+1)) newvort(i,j)=vort(i,j)+vincrement(i,j) temp=abs(vincrement(i,j)) if(temp.gt.vincmax)vincmax=temp 20 continue i=ix1 do 30 j=-iy1+1,channeltype*iy1-1 if(u(i,j).ge.0.)then au=-cx*u(i-1,j) bu= cx*u(i,j) cu= 0. else au=0 . bu=-cx*u(i,j) cu= cx*u(i-1,j) endif if(v(i,j).ge.0.)then av=-cy*v(i,j-1) bv= cy*v(i,j) cv= 0. else av= 0. bv=-cy*v(i,j) cv= cy*v(i,j+1) endif vincrement(i,j)= + -au*vort(i-1,j)-(bu+bv)*vort(i,j)-cu*vort(i-1,j) + -av*vort(i,j-1)-cv*vort(i,j+1) + +pex*(vort(i-1,j)-2*vort(i,j)+vort(i-1,j)) + +pey*(vort(i,j-1)-2*vort(i,j)+vort(i,j+1)) newvort(i,j)=vort(i,j)+vincrement(i,j) temp=abs(vincrement(i,j)) if(temp.gt.vincmax)vincmax=temp 30 continue return end c------------------------------------------------------------------ subroutine stepquick include 'step.h' vincmax=0. do 10 i=1,ixc1 do 10 j=1,iy1-1 call quick1(i,j,vinc) newvort(i,j)=vort(i,j)+vinc temp=abs(vinc) if(temp.gt.vincmax)vincmax=temp 10 continue do 20 i=ixc1+1,ix1-1 do 20 j=-iy1+1,channeltype*iy1-1 call quick2(i,j,vinc) newvort(i,j)=vort(i,j)+vinc temp=abs(vinc) if(temp.gt.vincmax)vincmax=temp 20 continue i=ix1 do 30 j=-iy1+1,channeltype*iy1-1 call quick3(i,j,vinc) newvort(i,j)=vort(i,j)+vinc temp=abs(vinc) if(temp.gt.vincmax)vincmax=temp 30 continue return end c---------------------------------------------------------------------- subroutine quick1(i,j,vortinc) include "step.h" c in inlet section, i<= ixc1 if(i.eq.1)then wleft=vort(0,j) else wleft=vort(i-2,j) endif if(j.eq.1)then wsouth=2*vort(i,j-1)-vort(i,j) else wsouth=vort(i,j-2) endif if(j.eq.(iy1-1))then wnorth=2*vort(i,j+1)-vort(i,j) else wnorth=vort(i,j+2) endif u0=u(i,j) v0=v(i,j) cox=u0*cx coy=v0*cy cox2=cox*cox coy2=coy*coy if(u0.ge.0.)then wxxx=(vort(i+1,j)-3*vort(i,j)+3*vort(i-1,j)-wleft)/6 else wxxx=(vort(i+2,j)-3*vort(i+1,j)+3*vort(i,j)-vort(i-1,j))/6 endif if(v0.ge.0.)then wyyy=(vort(i,j+1)-3*vort(i,j)+3*vort(i,j-1)-wsouth)/6 else wyyy=(wnorth-3*vort(i,j+1)+3*vort(i,j)-vort(i,j-1))/6 endif wxxy=((vort(i+1,j+1)-2*vort(i,j+1)+vort(i-1,j+1))- + (vort(i+1,j-1)-2*vort(i,j-1)+vort(i-1,j-1)))/4 wxyy=((vort(i+1,j+1)-2*vort(i+1,j)+vort(i+1,j-1))- + (vort(i-1,j+1)-2*vort(i-1,j)+vort(i-1,j-1)))/4 wxx=(vort(i+1,j)-2*vort(i,j)+vort(i-1,j))/2 wyy=(vort(i,j+1)-2*vort(i,j)+vort(i,j-1))/2 wxy=(vort(i+1,j+1)-vort(i+1,j-1)-vort(i-1,j+1)+vort(i-1,j-1))/4 wx=(vort(i+1,j)-vort(i-1,j))/2-wxxx wy=(vort(i,j+1)-vort(i,j-1))/2-wyyy vortinc=-cox*wx-coy*wy+(2*pex+cox2)*wxx+(2*pey+coy2)*wyy+ + cox*coy*wxy-(6*pex+cox2)*cox*wxxx-coy*(2*pex+cox2)*wxxy- + cox*(2*pey+coy2)*wxyy-(6*pey+coy2)*coy*wyyy return end subroutine quick2(i,j,vortinc) include "step.h" c in channel section, i> ixc1 wleft=vort(i-2,j) if(i.eq.(ixc1+1).and.((j.lt.0).or.(j.gt.iy1)))then wleft=2*vort(i-1,j)-vort(i,j) endif if(i.lt.(ix1-1))then wright=vort(i+2,j) else wright=vort(i,j) endif if(j.eq.(1-iy1))then wsouth=2*vort(i,j-1)-vort(i,j) else wsouth=vort(i,j-2) endif if(j.eq.(channeltype*iy1-1))then wnorth=2*vort(i,j+1)-vort(i,j) else wnorth=vort(i,j+2) endif u0=u(i,j) v0=v(i,j) cox=u0*cx coy=v0*cy cox2=cox*cox coy2=coy*coy if(u0.ge.0.)then wxxx=(vort(i+1,j)-3*vort(i,j)+3*vort(i-1,j)-wleft)/6 else wxxx=(wright-3*vort(i+1,j)+3*vort(i,j)-vort(i-1,j))/6 endif if(v0.ge.0.)then wyyy=(vort(i,j+1)-3*vort(i,j)+3*vort(i,j-1)-wsouth)/6 else wyyy=(wnorth-3*vort(i,j+1)+3*vort(i,j)-vort(i,j-1))/6 endif wxxy=((vort(i+1,j+1)-2*vort(i,j+1)+vort(i-1,j+1))- + (vort(i+1,j-1)-2*vort(i,j-1)+vort(i-1,j-1)))/4 wxyy=((vort(i+1,j+1)-2*vort(i+1,j)+vort(i+1,j-1))- + (vort(i-1,j+1)-2*vort(i-1,j)+vort(i-1,j-1)))/4 wxx=(vort(i+1,j)-2*vort(i,j)+vort(i-1,j))/2 wyy=(vort(i,j+1)-2*vort(i,j)+vort(i,j-1))/2 wxy=(vort(i+1,j+1)-vort(i+1,j-1)-vort(i-1,j+1)+vort(i-1,j-1))/4 wx=(vort(i+1,j)-vort(i-1,j))/2-wxxx wy=(vort(i,j+1)-vort(i,j-1))/2-wyyy vortinc=-cox*wx-coy*wy+(2*pex+cox2)*wxx+(2*pey+coy2)*wyy+ + cox*coy*wxy-(6*pex+cox2)*cox*wxxx-coy*(2*pex+cox2)*wxxy- + cox*(2*pey+coy2)*wxyy-(6*pey+coy2)*coy*wyyy return end subroutine quick3(i,j,vortinc) include "step.h" c in channel section, i=ix1 wleft=vort(i-2,j) wright=vort(i-2,j) if(j.eq.(1-iy1))then wsouth=2*vort(i,j-1)-vort(i,j) else wsouth=vort(i,j-2) endif if(j.eq.(channeltype*iy1-1))then wnorth=2*vort(i,j+1)-vort(i,j) else wnorth=vort(i,j+2) endif u0=u(i,j) v0=v(i,j) cox=u0*cx coy=v0*cy cox2=cox*cox coy2=coy*coy if(u0.ge.0.)then wxxx=(vort(i-1,j)-3*vort(i,j)+3*vort(i-1,j)-wleft)/6 else wxxx=(wright-3*vort(i-1,j)+3*vort(i,j)-vort(i-1,j))/6 endif if(v0.ge.0.)then wyyy=(vort(i,j+1)-3*vort(i,j)+3*vort(i,j-1)-wsouth)/6 else wyyy=(wnorth-3*vort(i,j+1)+3*vort(i,j)-vort(i,j-1))/6 endif wxxy=((vort(i-1,j+1)-2*vort(i,j+1)+vort(i-1,j+1))- + (vort(i-1,j-1)-2*vort(i,j-1)+vort(i-1,j-1)))/4 wxyy=((vort(i-1,j+1)-2*vort(i-1,j)+vort(i-1,j-1))- + (vort(i-1,j+1)-2*vort(i-1,j)+vort(i-1,j-1)))/4 wxx=(vort(i-1,j)-2*vort(i,j)+vort(i-1,j))/2 wyy=(vort(i,j+1)-2*vort(i,j)+vort(i,j-1))/2 wxy=(vort(i-1,j+1)-vort(i-1,j-1)-vort(i-1,j+1)+vort(i-1,j-1))/4 wx=(vort(i-1,j)-vort(i-1,j))/2-wxxx wy=(vort(i,j+1)-vort(i,j-1))/2-wyyy vortinc=-cox*wx-coy*wy+(2*pex+cox2)*wxx+(2*pey+coy2)*wyy+ + cox*coy*wxy-(6*pex+cox2)*cox*wxxx-coy*(2*pex+cox2)*wxxy- + cox*(2*pey+coy2)*wxyy-(6*pey+coy2)*coy*wyyy return end c------------------------------------------------------------------ subroutine calculatepsi include 'step.h' do 10 i=0,ix1 do 10 j=-iy1,channeltype*iy1 r1(i,j)=0.-newvort(i,j) 10 continue call mgrid(psi,r1,residual) return end c---------------------------------------------------------------------- subroutine calculateboundaryvorticity() include 'step.h' real *8 gxc,gyc gy=2./(dy*dy) gx=2./(dx*dx) gxc=2.*gx/3. gyc=2.*gy/3. do 20 i=1,ixc1-1 newvort(i,0)=gy*(psi(i,0)-psi(i,1)) newvort(i,iy1)=gy*(psi(i,iy1)-psi(i,iy1-1)) 20 continue i=ixc1 newvort(i,0)=gyc*(psi(i,0)-psi(i,1))+gxc*(psi(i,0)-psi(i+1,0)) if(channeltype.eq.1)then gyc=gy gxc=0. endif newvort(i,iy1)= + gyc*(psi(i,iy1)-psi(i,iy1-1))+gxc*(psi(i,iy1)-psi(i+1,iy1)) do 30 j=-iy1+1,-1 newvort(i,j)=gx*(psi(i,j)-psi(i+1,j)) 30 continue do 35 j=iy1+1,channeltype*iy1-1 newvort(i,j)=gx*(psi(i,j)-psi(i+1,j)) 35 continue do 40 i=ixc1+1,ix1 newvort(i,-iy1)=gy*(psi(i,-iy1)-psi(i,-iy1+1)) newvort(i,channeltype*iy1)= + gy*(psi(i,channeltype*iy1)-psi(i,channeltype*iy1-1)) 40 continue return end c------------------------------------------------------------------ subroutine calculatevelocities() include 'step.h' gx=1./dx gy=1./dy do 10 i=1,ixc1 do 10 j=1,iy1-1 u(i,j)=gy*(psi(i,j+1)-psi(i,j-1))/2 v(i,j)=gx*(psi(i-1,j)-psi(i+1,j))/2 10 continue do 20 i=ixc1+1,ix1-1 do 20 j=-iy1+1,channeltype*iy1-1 u(i,j)=gy*(psi(i,j+1)-psi(i,j-1))/2 v(i,j)=gx*(psi(i-1,j)-psi(i+1,j))/2 20 continue i=ix1 do 30 j=-iy1+1,channeltype*iy1-1 u(i,j)=gy*(psi(i,j+1)-psi(i,j-1))/2 v(i,j)=gx*(psi(i-1,j)-psi(i-1,j))/2 30 continue return end c-------------------------------------------------------------------- subroutine setflux() include 'step.h' if(currenttime.lt.0.25)then f=sin(2*3.14159*currenttime) else f=mean+(1.-mean)*sin(2*3.14159*currenttime) endif do 20 i=0,ixc1 psi(i,0)=0.-f psi(i,iy1)=f 20 continue do 30 j=-iy1,0 psi(ixc1,j)=0.-f 30 continue do 35 j=iy1,channeltype*iy1 psi(ixc1,j)=f 35 continue do 40 i=ixc1+1,ix1 psi(i,-iy1)=0.-f psi(i,channeltype*iy1)=f 40 continue do 50 j=0,iy1 y=j*dy u(0,j)=12*f*y*(1.-y) v(0,j)=0. newvort(0,j)=-12*f*(1-2*y) psi(0,j)=f*(6*y*y-4*y*y*y-1) 50 continue return end c-------------------------------------------------------------------- subroutine mgrid(z1,r1,difference) include 'step.h' difference=1. k=0 do while((difference.gt.tolerance).and.(k.lt.20)) call gauss (z1,r1,ixm1,iym1,iyp1,ix1,iy1,1,1) call residual(z1,r1,ixm1,iym1,iyp1,ix1,iy1, + r2,ixm2,iym2,iyp2,ix2,iy2,1 ) call zero (z2, ixm2,iym2,iyp2,ix2,iy2 ) call gauss (z2,r2,ixm2,iym2,iyp2,ix2,iy2,1,2) call residual(z2,r2,ixm2,iym2,iyp2,ix2,iy2, + r3,ixm3,iym3,iyp3,ix3,iy3,2 ) call zero (z3, ixm3,iym3,iyp3,ix3,iy3 ) if(iy1.gt.8)then call gauss (z3,r3,ixm3,iym3,iyp3,ix3,iy3,1,4) call residual(z3,r3,ixm3,iym3,iyp3,ix3,iy3, + r4,ixm4,iym4,iyp4,ix4,iy4,4 ) call zero (z4, ixm4,iym4,iyp4,ix4,iy4 ) call gauss (z4,r4,ixm4,iym4,iyp4,ix4,iy4,4,8) call interp (z3, ixm3,iym3,iyp3,ix3,iy3, + z4, ixm4,iym4,iyp4,ix4,iy4 ) endif call gauss (z3,r3,ixm3,iym3,iyp3,ix3,iy3,2,4) call interp (z2, ixm2,iym2,iyp2,ix2,iy2, + z3, ixm3,iym3,iyp3,ix3,iy3 ) call gauss (z2,r2,ixm2,iym2,iyp2,ix2,iy2,2,2) call interp (z1, ixm1,iym1,iyp1,ix1,iy1, + z2, ixm2,iym2,iyp2,ix2,iy2 ) call gauss (z1,r1,ixm1,iym1,iyp1,ix1,iy1,2,1) call diff (z1,r1,ixm1,iym1,iyp1,ix1,iy1,difference) k=k+1 end do if(k.eq.20)then write(screen,30) 30 format('warning: multigrid not converged') endif return end c--------------------------------------------------------------------- subroutine gauss(z,r,ixm,iym,iyp,ix,iy,repeats,level) implicit real*8(a-h,o-z) implicit integer (i-n) common/geometry/dx,dy,vfac,jv,ixc1,channeltype integer repeats,channeltype real *8 z(0:ixm,iym:iyp),r(0:ixm,iym:iyp) ixc=ixc1/level gx=1./(level*dx*level*dx) gy=1./(level*dy*level*dy) gxy=1./(2*(gx+gy)) do 30 k=1,repeats do 10 i=1,ixc do 10 j=1,iy-1 z(i,j)=gxy*( gx*(z(i+1,j)+z(i-1,j))+ + gy*(z(i,j+1)+z(i,j-1))-r(i,j) ) 10 continue do 20 i=ixc+1,ix-1 do 20 j=-iy+1,channeltype*iy-1 z(i,j)=gxy*( gx*(z(i+1,j)+z(i-1,j))+ + gy*(z(i,j+1)+z(i,j-1))-r(i,j) ) 20 continue i=ix do 25 j=-iy+1,channeltype*iy-1 z(i,j)=gxy*( gx*(z(i-1,j)+z(i-1,j))+ + gy*(z(i,j+1)+z(i,j-1))-r(i,j) ) 25 continue 30 continue return end c---------------------------------------------------------------------- subroutine zero(z,ixm,iym,iyp,ix,iy) implicit integer (i-n) real *8 z(0:ixm,iym:iyp) do 10 i=0,ix do 10 j=-iy,2*iy z(i,j)=0. 10 continue return end c---------------------------------------------------------------------- subroutine interp(z1,ixm1,iym1,iyp1,ix1,iy1, + z2,ixm2,iym2,iyp2,ix2,iy2) implicit integer (i-n) real *8 z1(0:ixm1,iym1:iyp1),z2(0:ixm2,iym2:iyp2) do 10 i=0,ix2 do 10 j=-iy2+1,2*iy2-1 z1(2*i,2*j)=z1(2*i,2*j)+z2(i,j) 10 continue do 20 i=1,ix1-1,2 do 20 j=-iy1,2*iy1,2 z1(i,j)=z1(i,j)+0.5*(z2((i+1)/2,j/2)+z2((i-1)/2,j/2)) 20 continue do 30 i=0,ix1,2 do 30 j=-iy1+1,2*iy1-1,2 z1(i,j)=z1(i,j)+0.5*(z2(i/2,(j+1)/2)+z2(i/2,(j-1)/2)) 30 continue do 40 i=1,ix1-1,2 do 40 j=-iy1+1,2*iy1-1,2 z1(i,j)=z1(i,j)+0.25*(z2((i+1)/2,(j+1)/2) + +z2((i+1)/2,(j-1)/2)+z2((i-1)/2,(j+1)/2) + +z2((i-1)/2,(j-1)/2)) 40 continue return end c----------------------------------------------------------------------- subroutine residual(z1,r1,ixm1,iym1,iyp1,ix1,iy1, + r2,ixm2,iym2,iyp2,ix2,iy2,level) implicit real *8(a-h,o-z) implicit integer (i-n) common/geometry/dx,dy,vfac,jv,ixc1,channeltype integer channeltype real *8 z1(0:ixm1,iym1:iyp1),r1(0:ixm1,iym1:iyp1), + r2(0:ixm2,iym2:iyp2) gx=1./(level*dx*level*dx) gy=1./(level*dy*level*dy) do 20 j=-iy1+2,2*iy1-2,2 do 10 i=2,ix1-2,2 r2(i/2,j/2)=r1(i,j)-gx*(z1(i+1,j)-2*z1(i,j)+z1(i-1,j)) + -gy*(z1(i,j+1)-2*z1(i,j)+z1(i,j-1)) 10 continue i=ix1 r2(i/2,j/2)=r1(i,j)-gx*(z1(i-1,j)-2*z1(i,j)+z1(i-1,j)) + -gy*(z1(i,j+1)-2*z1(i,j)+z1(i,j-1)) 20 continue return end c----------------------------------------------------------------------- subroutine diff(z1,r1,ixm1,iym1,iyp1,ix1,iy1,w) implicit real *8(a-h,o-z) implicit integer (i-n) common/geometry/dx,dy,vfac,jv,ixc1,channeltype integer channeltype real *8 z1(0:ixm1,iym1:iyp1),r1(0:ixm1,iym1:iyp1) gx=1./(dx*dx) gy=1./(dy*dy) w=0. d=0. do 10 i=1,ixc1 do 10 j=1,iy1-1 d= r1(i,j)- + gx*( z1(i+1,j)-2*z1(i,j)+z1(i-1,j) ) + -gy*( z1(i,j+1)-2*z1(i,j)+z1(i,j-1) ) w=w+abs(d) 10 continue do 20 i=ixc1+1,ix1-1 do 20 j=-iy1+1,channeltype*iy1-1 d= r1(i,j)- + gx*( z1(i+1,j)-2*z1(i,j)+z1(i-1,j) ) + -gy*( z1(i,j+1)-2*z1(i,j)+z1(i,j-1) ) w=w+abs(d) 20 continue i=ix1 do 30 j=-iy1+1,channeltype*iy1-1 d= r1(i,j)- + gx*( z1(i-1,j)-2*z1(i,j)+z1(i-1,j) ) + -gy*( z1(i,j+1)-2*z1(i,j)+z1(i,j-1) ) w=w+abs(d) 30 continue d=abs(z1(0,iy1)-z1(0,0)) if(d.gt.0.)then w=w/(d*ix1*iy1) else w=w/(ix1*iy1) endif return end c---------------------------------------------------------------------- subroutine start() include 'step.h' keyboard =5 screen =6 arcfile =8 logfile =9 inputfile=10 sumfile =11 psifile =12 vortfile =13 open(unit=inputfile,file='chan.dat',form='formatted') open(unit=logfile ,file='chan.log',form='formatted') open(unit=psifile ,file='chan.psi',form='formatted') open(unit=vortfile ,file='chan.vor',form='formatted') open(unit=sumfile ,file='chan.sum',form='formatted') c if method>4 then start from archive with method-5 as method central =0 upwind =1 ncupwind =2 laxwendroff =3 quickest =4 finish =0 swapflow =1 changereynolds =2 changedt =3 changemethod =4 archiveflow =5 read(inputfile,1)channeltype read(inputfile,1)method read(inputfile,1)ix1 read(inputfile,1)iy1 read(inputfile,1)ixc1 read(inputfile,1)ns1 read(inputfile,1)ns2 read(inputfile,2)reynolds read(inputfile,2)strouhal read(inputfile,2)tolerance read(inputfile,2)mean read(inputfile,2)dt 1 format(i5) 2 format(f10.3) ix2=ix1/2 ix3=ix2/2 ix4=ix3/2 iy2=iy1/2 iy3=iy2/2 iy4=iy3/2 dy=1./iy1 dx=dy jv=6.38/dx vfac=(6.38-jv*dx)/dx jv=jv+ixc1 cx=dt/(strouhal*dx) cy=dt/(strouhal*dy) pex=dt/(reynolds*strouhal*dx*dx) pey=dt/(reynolds*strouhal*dy*dy) if(channeltype.eq.1)then write(screen, 3) write(logfile,3) elseif(channeltype.eq.2)then write(screen,4) write(logfile,4) else write(screen,5) write(logfile,5) endif 3 format('Backstep') 4 format('Symmetric Expansion') 5 format('Panic: unknown channel type') if((method.eq.0).or.((method-5).eq.0))then write(screen,6) write(logfile,6) else if((method.eq.1).or.((method-5).eq.1))then write(screen,7) write(logfile,7) else if((method.eq.2).or.((method-5).eq.2))then write(screen,8) write(logfile,8) else if((method.eq.3).or.((method-5).eq.3))then write(screen,9) write(logfile,9) else if((method.eq.4).or.((method-5).eq.4))then write(screen,10) write(logfile,10) else write(screen,11) write(logfile,11) endif 6 format('Method: central difference') 7 format('Method: upwind difference') 8 format('Method: nonconservative upwind difference') 9 format('Method: Lax Wendroff difference') 10 format('Method: Quickest difference') 11 format('Method: Panic: unkown method') write(screen,20)ns1,ns2,reynolds,strouhal,tolerance,mean,dt write(logfile,20)ns1,ns2,reynolds,strouhal,tolerance,mean,dt 20 format(i6,'steps then',i5,' steps per print',/, + 'reynolds number = ',f8.1,/, + 'strouhal number = ',f8.4,/, + 'tolerance = ',f8.6,/, + 'mean flow = ',f8.3,/, + 'time step = ',f8.6) write(screen,30)dx*ix1,dx*(ix1-ixc1),cx,cy,pex,pey write(logfile,30)dx*ix1,dx*(ix1-ixc1),cx,cy,pex,pey 30 format('channel length',f10.2,/, + 'lee length ',f10.2,/, + 'cx,cy ',2f10.4,/, + 'pex,pey ',2f10.4) if(method.ge.5)then method=method-5 call unarchive() else currenttime=0. call initialiseall endif return end c--------------------------------------------------------------------- c---------------------------------------------------------------------- subroutine archive() include 'step.h' open(unit=arcfile ,file='chan.arc',form='unformatted') write(arcfile)channeltype write(arcfile)method write(arcfile)ix1 write(arcfile)iy1 write(arcfile)ixc1 write(arcfile)ns1 write(arcfile)ns2 write(arcfile)reynolds write(arcfile)strouhal write(arcfile)tolerance write(arcfile)mean write(arcfile)dt write(arcfile)currenttime write(arcfile)psi write(arcfile)oldvort write(arcfile)vort close(unit=arcfile) return end c---------------------------------------------------------------------- subroutine unarchive() include 'step.h' integer channeltypearc,methodarc,ix1arc,iy1arc,ixc1arc, + ns1arc,ns2arc real*8 reynoldsarc,strouhalarc,tolerancearc,dtarc c c want to be able to unarchive from coarser resolutions c but assuming global array sizes are same. Use oldvort c as dummy array in this case c open(unit=arcfile ,file='chan.arc',form='unformatted') read(arcfile)channeltypearc read(arcfile)methodarc read(arcfile)ix1arc read(arcfile)iy1arc read(arcfile)ixc1arc read(arcfile)ns1arc read(arcfile)ns2arc read(arcfile)reynoldsarc read(arcfile)strouhalarc read(arcfile)tolerancearc read(arcfile)meanarc read(arcfile)dtarc read(arcfile)currenttime if(iy1arc.eq.iy1)then write(screen,5) 5 format('unarchiving normal data') read(arcfile)psi read(arcfile)oldvort read(arcfile)vort else if(2*iy1arc.eq.iy1)then write(screen,8) 8 format('unarchiving from coarse data') read(arcfile)oldvort c 1)set values from coarse array do 10 i=0,ix1arc do 10 j=-iy1arc,2*iy1arc psi(2*i,2*j)=oldvort(i,j) 10 continue c 2)interpolate cell by cell do 20 i=0,ix1-2,2 do 20 j=-iy1,2*iy1-2,2 psi(i+1,j)=0.5*(psi(i,j)+psi(i+2,j)) psi(i,j+1)=0.5*(psi(i,j)+psi(i,j+2)) psi(i+1,j+1)=0.25*(psi(i,j)+psi(i+2,j)+ + psi(i,j+2)+psi(i+2,j+2)) 20 continue do 30 i=0,ix1-2,2 psi(i+1,2*iy1)=0.5*(psi(i,2*iy1)+psi(i+2,2*iy1)) 30 continue do 40 j=-iy1,2*iy1-2,2 psi(ix1,j+1)=0.5*(psi(ix1,j)+psi(ix1,j+2)) 40 continue read(arcfile)oldvort read(arcfile)oldvort c 1)set values from coarse array do 50 i=0,ix1arc do 50 j=-iy1arc,2*iy1arc vort(2*i,2*j)=oldvort(i,j) 50 continue c 2)interpolate cell by cell do 60 i=0,ix1-2,2 do 60 j=-iy1,2*iy1-2,2 vort(i+1,j)=0.5*(vort(i,j)+vort(i+2,j)) vort(i,j+1)=0.5*(vort(i,j)+vort(i,j+2)) vort(i+1,j+1)=0.25*(vort(i,j)+vort(i+2,j)+ + vort(i,j+2)+vort(i+2,j+2)) 60 continue do 70 i=0,ix1-2,2 vort(i+1,2*iy1)=0.5*(vort(i,2*iy1)+vort(i+2,2*iy1)) 70 continue do 80 j=-iy1,2*iy1-2,2 vort(ix1,j+1)=0.5*(vort(ix1,j)+vort(ix1,j+2)) 80 continue do 90 i=0,ix1 do 90 j=-iy1,2*iy1 oldvort(i,j)=vort(i,j) 90 continue else if(2*iy1.eq.iy1arc)then write(screen,92) 92 format('unarchiving fine data') read(arcfile)oldvort do 100 i=0,ix1 do 100 j=-iy1,2*iy1 psi(i,j)=oldvort(2*i,2*j) 100 continue read(arcfile)newvort do 110 i=0,ix1 do 110 j=-iy1,2*iy1 oldvort(i,j)=newvort(2*i,2*j) 110 continue read(arcfile)newvort do 120 i=0,ix1 do 120 j=-iy1,2*iy1 vort(i,j)=newvort(2*i,2*j) 120 continue else write(screen,95)ix1,iy1,ix1arc,iy1arc 95 format('PANIC in restorearc:',4i6) endif call calculatevelocities close(unit=arcfile) return end