program square c********************************************************************** c c MAIN PROGRAM c c solve the navier-stokes equations for two incompressible fluids c in a rectangular region with a uniform gird. The two fluids are c separated by an interface c c parameters: c c tend: end of calculation c rey: Reynolds number c cfl: cfl number c (nx,ny): grid size c c c********************************************************************** integer nx,ny,nxy,karea,nxm,nym,nm,nme,mxy,mxym, + nrew,npit,kpit,nout,i,nstep,kk,kunit,krei logical piter include 'gridsz' parameter (nrew=20,npit=2,karea=5,krei=10) parameter (nm=nxm*nym, nme=(nxm+1)*(nym+1)) double precision uf(nm),vf(nm),px(nm),py(nm), + pf(0:nxm+1,0:nym+1) double precision us(nm),vs(nm), + fx(nm),fy(nm),wp(nm),wy(nm), + ux(nme),uy(nme),vx(nme),vy(nme) double precision tout(10) double precision pi,tend,rey,cfl,h,hinv,hhinv,dtin, + dt,t,ek,dth,ab,rin,ru,rb,grav common/reynod/rey common/meshsz/h,hinv common/hfract/hhinv common/timest/dt common/timspc/dth common/valupi/pi common/gravity/grav data kpit,kunit/0,20/, t,ru,rb,grav/0.D0,1.D0,0.5D0,10.D0/ data piter/.true./ c ********************************************************* read (5,*) nx,ny read (5,*) tend,rey,cfl read (5,*) nout read (5,*) (tout(i),i=1,nout) pi = 4.D0*datan(1.D0) nxy = nx*ny mxy = 2*nxy hinv = dble(ny) hhinv = 0.5D0*hinv h = 1.D0/hinv dtin = 0.5D0*h nstep = idint((tend+1.D-10)/dtin) rin = 0.1D0 c initialize the flow call stflow(uf,vf,px,py,nx,ny) call inifrt(pf,rin,nx,ny) call densit(pf,ru,rb,nx,ny) call output(pf,uf,vf,nx,ny,kunit) c******************************************************************* c c--- -- time evolution -- c c do 5 kk=1,nstep+npit kk = 0 5 continue call chtmst(uf,vf,ek,dtin,cfl,nx,ny) dth = dt*hinv ab = 0.5D0*dt*hinv*hinv/rey call edgevl(uf,vf,us,vs,ux,vx,uy,vy,nx,ny) call macpro(ux,vx,uy,vy,nx,ny,mxy) call timrhs(uf,vf,px,py,us,vs,ux,vx,uy,vy,fx,fy,ab,nx,ny) if( .not. piter ) call convec(uf,vf,ux,vy,pf,nx,ny) if( mod(kk,krei) .eq. 0 ) call reinit(pf,nx,ny) c--- obtain the intermediate velocity u star call conjug(fx,us,uf,wp,wy,ab,nxy,nx,ny) call conjug(fy,vs,vf,wp,wy,ab,nxy,nx,ny) c--- project to the incompressible velocity space call projec(uf,vf,us,vs,px,py,nx,ny,mxy,piter) c if( piter ) then kpit = kpit + 1 if( kpit .eq. npit) piter = .false. else t = t + dt kk = kk + 1 end if print 110, dt,t,ek c ......................................................... do 150 i=1,nout if ( t-tout(i) .ge. 0. .and. t-tout(i) .lt. dt) then kunit = kunit+1 c call vortic(uf,vf,wp,nx,ny) call output(pf,uf,vf,nx,ny,kunit) c call stream(uf,vf,vx,nx,ny) endif 150 continue if ( t .le. tend) goto 5 c 5 continue c ************************************************* 110 format(' dt=',f10.7,' t=',f8.5,' ek=',f12.8 ) end subroutine densit(pf,ru,rb,nx,ny) integer i,j,nxm,nym,mxym,nx,ny include 'gridsz' double precision pf(0:nx+1,0:ny+1) double precision roa,ru,rb common/density/roa(0:nxm+1,0:nym+1) do 10 j=0,ny+1 do 10 i=0,nx+1 if (pf(i,j) .ge. 0.D0) then roa(i,j) = rb else roa(i,j) = ru endif 10 continue end subroutine stflow(uf,vf,px,py,nx,ny) *********************************************************** c c initialize the flow by specifying the velocity and c set pressure to zero. c c uf,vf <- initial value to be set c c********************************************************** integer i,j,nx,ny,nxm,nym,mxym include 'gridsz' double precision uf(nx,ny),vf(nx,ny),px(nx,ny),py(nx,ny) double precision h,hinv,x,y,pi,uby,vby,utop,ubot common/meshsz/h,hinv common/uvboud/uby(nym),vby(nym),utop,ubot c ****************************************************** pi = 4.D0*datan(1.D0) do 5 j=1,ny y = (dble(j)-0.5D0)*h do 5 i=1,nx x = (dble(i)-0.5D0)*h uf(i,j) = 0.D0 vf(i,j) = 0.D0 px(i,j) = 0.D0 py(i,j) = 0.D0 5 continue do 6 j=1,ny uby(j) = 0.D0 vby(j) = 0.D0 6 continue utop = 0.D0 ubot = 0.D0 * ***************************************************** end subroutine stream(uf,vf,psi,nx,ny) integer i,j,nx,ny,nxm,nym,mxym include 'gridsz' double precision uf(nx,ny),vf(nx,ny),psi(0:nx,0:ny) double precision uby,vby,utop,ubot,h,hinv common/uvboud/uby(nym),vby(nym),utop,ubot common/meshsz/h,hinv do 5 i=0,nx psi(i,0) = 0.D0 5 continue do 10 j=1,ny psi(0,j) = 0.D0 do 12 i=1,nx psi(i,j) = psi(i,j-1) - psi(i-1,j) + psi(i-1,j-1) + + 2.D0*h*uf(i,j) 12 continue 10 continue write(24,*) psi write(24,101) 101 format(/) end subroutine edgevl(u,v,ud,vd,ux,vx,uy,vy,nx,ny) c**************************************************************** c c determine the edge values of velocity by upwinding, c before making the pressure correction. c c (u,v) -> velocity field c (ux,vx) (uy,vy) <- 1), constructed first derivatives c <- 2), edge values c (ud,vd) <- laplacian of (u,v) c temp space: c (up,vp) (uq,vq) <- first derivatives after c 2nd order upwinding used c in transverse directions. c c**************************************************************** integer i,j,nx,ny,nxm,nym,mxym include 'gridsz' double precision u(nx,ny),v(nx,ny), + ud(nx,ny),vd(nx,ny), + ux(0:nx,0:ny),vx(0:nx,0:ny), + uy(0:nx,0:ny),vy(0:nx,0:ny) double precision up(nxm,nym),vp(nxm,nym), + uq(nxm,nym),vq(nxm,nym) double precision uld,vld,h,hinv,rey,ul,ur,vl,vr,ub,ut,vb,vt, + roa,hf,hdt,dt,uby,vby,utop,ubot,grav parameter (hf=0.5D0) common/timest/dt common/meshsz/h,hinv common/reynod/rey common/density/roa(0:nxm+1,0:nym+1) common/gravity/grav common/uvboud/uby(nym),vby(nym),utop,ubot c initialize the local arrays: do 2 j=1,ny do 2 i=1,nx up(i,j) = 0.D0 vp(i,j) = 0.D0 uq(i,j) = 0.D0 vq(i,j) = 0.D0 2 continue c hdt = hf*dt c c approximate the first derivatives and laplacian of velocities: call evalim(u,ux,uy,uby,utop,ubot,nx,ny) call evalim(v,vx,vy,vby,0.D0,0.D0,nx,ny) call evalap(u,uby,utop,ubot,ud,nx,ny) call evalap(v,vby,0.D0,0.D0,vd,nx,ny) c c *************************************************************** c construct the tangential derivatives along cell edges c c 1): along vertical edges, do 10 j=1,ny up(1,j) = ux(1,j) vp(1,j) = vx(1,j) up(nx,j) = ux(nx,j) vp(nx,j) = vx(nx,j) 10 continue do 20 j=1,ny do 20 i=2,nx-1 call upwsec( u(i,j), u(i-1,j), u(i+1,j), + v(i,j), v(i-1,j), v(i+1,j), + ux(i,j), ux(i-1,j), ux(i+1,j), + vx(i,j), vx(i-1,j), vx(i+1,j), + up(i,j), vp(i,j) ) 20 continue c 2): along horizontal edges, do 40 i=1,nx uq(i,1) = uy(i,1) vq(i,1) = vy(i,1) uq(i,ny) = uy(i,ny) vq(i,ny) = vy(i,ny) 40 continue do 60 j=2,ny-1 do 60 i=1,nx call upwsec( v(i,j), v(i,j-1), v(i,j+1), + u(i,j), u(i,j-1), u(i,j+1), + vy(i,j), vy(i,j-1), vy(i,j+1), + uy(i,j), uy(i,j-1), uy(i,j+1), + vq(i,j), uq(i,j) ) 60 continue c *********************************************************** c construct time derivatives c do 70 j=1,ny do 70 i=1,nx uld = ud(i,j)/rey/roa(i,j) vld = vd(i,j)/rey/roa(i,j) uq(i,j) = uld - hinv*( u(i,j)*ux(i,j) + v(i,j)*uq(i,j) ) up(i,j) = uld - hinv*( u(i,j)*up(i,j) + v(i,j)*uy(i,j) ) vq(i,j) = vld - hinv*( u(i,j)*vx(i,j) + v(i,j)*vq(i,j) ) - grav vp(i,j) = vld - hinv*( u(i,j)*vp(i,j) + v(i,j)*vy(i,j) ) - grav 70 continue c *********************************************************** c Taylor extrapolation and upwinding to obtain edge values, c vertical: do 80 j=1,ny do 80 i=1,nx-1 ul = u(i,j) + hdt*uq(i,j) + ux(i,j)*hf ur = u(i+1,j) + hdt*uq(i+1,j) - ux(i+1,j)*hf vl = v(i,j) + hdt*vq(i,j) + vx(i,j)*hf vr = v(i+1,j) + hdt*vq(i+1,j) - vx(i+1,j)*hf call upwind(ul,ur,vl,vr,ux(i,j),vx(i,j)) 80 continue c horizontal: do 85 j=1,ny-1 do 85 i=1,nx ub = u(i,j) + hdt*up(i,j) + uy(i,j)*hf ut = u(i,j+1) + hdt*up(i,j+1) - uy(i,j+1)*hf vb = v(i,j) + hdt*vp(i,j) + vy(i,j)*hf vt = v(i,j+1) + hdt*vp(i,j+1) - vy(i,j+1)*hf call upwind(vb,vt,ub,ut,vy(i,j),uy(i,j)) 85 continue c--- expanding the array (uy,vy) and (ux,vx) by incorporating c boundary values. do 90 i=1,nx uy(i,0) = ubot vy(i,0) = 0.D0 uy(i,ny) = utop vy(i,ny) = 0.D0 90 continue do 95 j=1,ny ux(0,j) = uby(j) vx(0,j) = 0.D0 ux(nx,j) = uby(j) vx(nx,j) = 0.D0 95 continue end subroutine evalap(u,uby,utop,ubot,ud,nx,ny) c********************************************************************* c c compute the Laplacian of u and add its effect to the velocity. c central differencing in the interior and one-sided differencing c on the boundaries c c u -> u at cell centers c uby -> u at the vertical walls c ud <- Laplacian of u. c c********************************************************************* integer nx,ny,i,j double precision u(nx,ny),ud(nx,ny),uby(ny),utop,ubot, + c2,c3,cb,sv,tn,h,hinv parameter (c2=2.D0,c3=-0.2D0,cb=3.2D0,sv=7.D0,tn=10.D0) common/meshsz/h,hinv c .......................................................... c--- interior do 10 j=2,ny-1 do 10 i=2,nx-1 ud(i,j) = u(i+1,j) + u(i-1,j) + + u(i,j+1) + u(i,j-1) - 4.D0*u(i,j) 10 continue c .......................................................... c boundaries do 20 i=2,nx-1 ud(i,1) = u(i+1,1) + u(i-1,1) + cb*ubot + + c2*u(i,2) + c3*u(i,3) - sv*u(i,1) ud(i,ny) = u(i+1,ny) + u(i-1,ny) + cb*utop + + c2*u(i,ny-1) + c3*u(i,ny-2) - sv*u(i,ny) 20 continue do 30 j=2,ny-1 ud(1,j) = u(1,j+1) + u(1,j-1) + cb*uby(j) + + c2*u(2,j) + c3*u(3,j) - sv*u(1,j) ud(nx,j) = u(nx,j+1) + u(nx,j-1) + cb*uby(j) + + c2*u(nx-1,j) + c3*u(nx-2,j) - sv*u(nx,j) 30 continue c corners ud(1,1) = c2*(u(1,2)+u(2,1)) + + cb*(uby(1)+ubot) + + c3*(u(1,3)+u(3,1)) - + tn*u(1,1) ud(nx,1) = c2*(u(nx-1,1)+u(nx,2)) + + cb*(uby(1)+ubot) + + c3*(u(nx-2,1)+u(nx,3)) - + tn*u(nx,1) ud(1,ny) = c2*(u(1,ny-1)+u(2,ny)) + + cb*(uby(ny)+utop) + + c3*(u(1,ny-2)+u(3,ny)) - + tn*u(1,ny) ud(nx,ny) = c2*(u(nx-1,ny)+u(nx,ny-1)) + + cb*(uby(ny)+utop) + + c3*(u(nx-2,ny)+u(nx,ny-2)) - + tn*u(nx,ny) do 40 j=1,ny do 40 i=1,nx ud(i,j) = hinv*hinv*ud(i,j) 40 continue end subroutine evalim(u,ux,uy,uby,utop,ubot,nx,ny) c******************************************************************* c c construct the first derivatives with limiters c c u -> u at cell centers c ux,uy <- constructed first derivatives c c******************************************************************* integer i,j,nx,ny double precision u(nx,ny),ux(0:nx,0:ny),uy(0:nx,0:ny),uby(ny) double precision blimit,slimit,utop,ubot c .......................................................... c determine ux. do 10 j=1,ny ux(1,j) = blimit( uby(j), u(1,j), u(2,j) ) ux(nx,j) =-blimit( uby(j), u(nx,j), u(nx-1,j) ) 10 continue do 20 j=1,ny do 20 i=2,nx-1 ux(i,j) = slimit( u(i+1,j)-u(i,j), u(i,j)-u(i-1,j) ) 20 continue c .......................................................... c determine uy. do 60 i=1,nx uy(i,1) = blimit( ubot, u(i,1), u(i,2) ) uy(i,ny) =-blimit( utop, u(i,ny), u(i,ny-1) ) 60 continue do 70 j=2,ny-1 do 70 i=1,nx uy(i,j) = slimit( u(i,j+1)-u(i,j), u(i,j)-u(i,j-1) ) 70 continue end subroutine timrhs(uf,vf,px,py,ud,vd,ux,vx,uy,vy,fx,fy,ab,nx,ny) c******************************************************************* c c evaluate the forcing terms for the time stepping c c u + dt*( 1/2/rey grad^2 u - (u grad u) - grad p ) c c c******************************************************************* integer i,j,nx,ny,nxm,nym,mxym include 'gridsz' double precision fx(nx,ny),fy(nx,ny), + uf(nx,ny),vf(nx,ny), + ud(nx,ny),vd(nx,ny), + px(nx,ny),py(nx,ny), + ux(0:nx,0:ny),vx(0:nx,0:ny), + uy(0:nx,0:ny),vy(0:nx,0:ny), + dt,dth,hdtrey,rey,hf,ufx,vfx,uav,vav, + ab,uby,vby,utop,ubot,fadb,fadt,roa,grav common/timest/dt common/timspc/dth common/reynod/rey common/uvboud/uby(nym),vby(nym),utop,ubot common/gravity/grav common/density/roa(0:nxm+1,0:nym+1) parameter (hf=0.5D0) hdtrey = hf*dt/rey c ........................................................ do 80 j=1,ny do 80 i=1,nx uav = hf*( ux(i-1,j) + ux(i,j) ) vav = hf*( vy(i,j-1) + vy(i,j) ) ufx = uav*( ux(i,j) - ux(i-1,j) ) + + vav*( uy(i,j) - uy(i,j-1) ) vfx = uav*( vx(i,j) - vx(i-1,j) ) + + vav*( vy(i,j) - vy(i,j-1) ) fx(i,j) = uf(i,j) + (hdtrey*ud(i,j) - + dt*px(i,j))/roa(i,j) - dth*ufx fy(i,j) = vf(i,j) + (hdtrey*vd(i,j) - + dt*py(i,j))/roa(i,j) - dth*vfx - + dt*grav 80 continue do 90 j=1,ny fx(1,j) = fx(1,j) + 3.2D0*ab*uby(j) fx(nx,j) = fx(nx,j) + 3.2D0*ab*uby(j) 90 continue fadb = 3.2D0*ab*ubot fadt = 3.2D0*ab*utop do 95 i=1,nx fx(i,1) = fx(i,1) + fadb fx(i,ny) = fx(i,ny) + fadt 95 continue end subroutine conjug(rem,x,x0,p,y,ab,nxy,nx,ny) c********************************************************************* c c solve the linear system A x = rem c by the unconditioned conjugate gradient method c c p,y : work space c nxy -> dimension of rem and x c ab -> lambda in CN scheme matrix c c matrix multiplication by routine atimex c accuracy controlled by: c || e || 2 <= eps. c c********************************************************************* integer nxy,nx,ny,i,its,itmax double precision rem(nxy),x(nxy),x0(nxy),p(nxy),y(nxy), + ab,eps,rnosq,rnnsq,papt,ak,bk parameter (eps=1.D-10) c--- initialize do 5 i=1,nxy x(i) = x0(i) 5 continue call atimex(x,y,ab,nxy,nx,ny) itmax = nxy rnosq = 0.D+00 do 10 i=1,nxy rem(i) = rem(i)-y(i) p(i) = rem(i) rnosq = rnosq + rem(i)*rem(i) 10 continue c--- starting iterations do 50 its=1,itmax call atimex(p,y,ab,nxy,nx,ny) papt = 0.D+00 do 20 i=1,nxy papt = papt + p(i)*y(i) 20 continue if( papt .eq. 0.D+00 ) return ak = rnosq/papt rnnsq = 0.D+00 do 30 i=1,nxy x(i) = x(i)+ak*p(i) rem(i) = rem(i)-ak*y(i) rnnsq = rnnsq + rem(i)*rem(i) 30 continue if( dsqrt(rnnsq/dble(nxy)) .le. eps ) return bk = rnnsq/rnosq do 40 i=1,nxy p(i) = rem(i)+bk*p(i) 40 continue 50 rnosq=rnnsq pause 'maximum iterations exceeded' c write(6,*) 'CNCG =', its, ' rnnsq=',rnnsq end subroutine atimex(x,y,ab,nxy,nx,ny) c********************************************************************** c c matrix multiplication y = A x c c x -> vector of length nxy c ab -> 0.5*dt/h/h/rey c y <- vector of length nxy c c A = discretization of the operator (I-dt/2/rey*grad^2) c c********************************************************************** integer nxy,nx,ny,nxm,nym,mxym include 'gridsz' double precision x(nxy),y(nxy),ab,roa common/density/roa(0:nxm+1,0:nym+1) integer i,j,k,k1,k2 double precision tn,sv,c2,c3 parameter (sv=7.D0, tn=10.D0, c2=2.D0, c3=-0.2D0) c********************************************************** c interior points. k=nx+2 do 10 j=2,ny-1 do 15 i=2,nx-1 y(k) = x(k) + ab*(4.D0*x(k) - x(k-1) - x(k+1) - + x(k-nx) - x(k+nx) )/roa(i,j) k = k+1 15 continue k = k+2 10 continue c********************************************************** c boundary points. k1=2 k2=nx*(ny-1)+2 do 20 i=2,nx-1 y(k1) = x(k1) + ab*(sv*x(k1) - x(k1-1) - x(k1+1) - + c2*x(k1+nx) - c3*x(k1+2*nx) )/ + roa(i,1) y(k2) = x(k2) + ab*(sv*x(k2) - x(k2-1) - x(k2+1) - + c2*x(k2-nx) - c3*x(k2-2*nx) )/ + roa(i,ny) k1=k1+1 k2=k2+1 20 continue k1=nx+1 k2=2*nx do 30 j=2,ny-1 y(k1) = x(k1) + ab*(sv*x(k1) - x(k1-nx) - x(k1+nx) - + c2*x(k1+1) - c3*x(k1+2))/ + roa(1,j) y(k2) = x(k2) + ab*(sv*x(k2) - x(k2-nx) - x(k2+nx) - + c2*x(k2-1) - c3*x(k2-2))/ + roa(nx,j) k1=k1+nx k2=k2+nx 30 continue c*********************************************************** c corners k=1 y(k) = x(k) + ab*(tn*x(k) - c2*(x(k+1)+x(k+nx)) - + c3*(x(k+2)+x(k+2*nx)) )/roa(1,1) k=nx y(k) = x(k) + ab*(tn*x(k) - c2*(x(k-1)+x(k+nx)) - + c3*(x(k-2)+x(k+2*nx)))/roa(nx,1) k=nx*(ny-1)+1 y(k) = x(k) + ab*(tn*x(k) - c2*(x(k-nx)+x(k+1)) - + c3*(x(k-2*nx)+x(k+2)))/roa(1,ny) k=nx*ny y(k) = x(k) + ab*(tn*x(k) - c2*(x(k-1)+x(k-nx)) - + c3*(x(k-2)+x(k-2*nx)))/roa(nx,ny) end subroutine convec(uf,vf,uv,vh,pf,nx,ny) * ********************************************************** * * nonsplit second-order upwind scheme for the linear advection * in the backfacing step channel: * * phi_t + u phi_x + v phi_y = 0 * * velocity field given at cell centres and cell edges. * * ********************************************************** integer i,j,nx,ny,nxm,nym,mxym include 'gridsz' double precision dth,slimit double precision uf(nx,ny),vf(nx,ny), + uv(0:nx,0:ny),vh(0:nx,0:ny), + pf(0:nx+1,0:ny+1) double precision px(nxm,nym),py(nxm,nym), + pv(0:nxm,0:nym),ph(0:nxm,0:nym) common/timspc/dth * construct 1st derivatives with limiter in the x (normal) direction * and upwind in the y (transverse) direction do 5 j=1,ny do 5 i=1,nx px(i,j) = slimit(pf(i,j)-pf(i-1,j),pf(i+1,j)-pf(i,j)) if (vf(i,j) .gt. 0.D0) then py(i,j) = pf(i,j) - pf(i,j-1) else py(i,j) = pf(i,j+1) - pf(i,j) end if 5 continue * approximate vertical edge values. do 10 j=1,ny do 10 i=1,nx-1 if (uv(i,j) .gt. 0.D0) then pv(i,j) = pf(i,j) + 0.5D0*( (1.D0 - + dth*uf(i,j) )*px(i,j) - + dth*vf(i,j) *py(i,j) ) else pv(i,j) = pf(i+1,j) - 0.5D0*( (1.D0 + + dth*uf(i+1,j) )*px(i+1,j) + + dth*vf(i+1,j) *py(i+1,j) ) end if 10 continue do 22 j=1,ny pv(0,j) = pv(1,j) pv(nx,j) = pv(nx-1,j) 22 continue * * construct 1st derivatives with limiter in the y (normal) direction * and upwind in the x (transverse) direction do 25 j=1,ny do 25 i=1,ny py(i,j) = slimit(pf(i,j)-pf(i,j-1),pf(i,j+1)-pf(i,j)) if (uf(i,j) .gt. 0.D0) then px(i,j) = pf(i,j) - pf(i-1,j) else px(i,j) = pf(i+1,j) - pf(i,j) end if 25 continue do 30 j=1,ny-1 do 30 i=1,nx if (vh(i,j) .gt. 0.D0) then ph(i,j) = pf(i,j) + 0.5D+00*( (1.0D+00 - + dth*vf(i,j) )*py(i,j) - + dth*uf(i,j) *px(i,j) ) else ph(i,j) = pf(i,j+1) - 0.5D+00*( (1.0D+00 + + dth*vf(i,j+1) )*py(i,j+1) + + dth*uf(i,j+1) *px(i,j+1) ) end if 30 continue do 42 i=1,nx ph(i,0) = ph(i,1) ph(i,ny) = ph(i,ny-1) 42 continue * * calculate fluxes and update pf. do 50 j=1,ny do 50 i=1,nx pf(i,j) = pf(i,j) + dth*( + uv(i-1,j)*pv(i-1,j) - uv(i,j)*pv(i,j) + + vh(i,j-1)*ph(i,j-1) - vh(i,j)*ph(i,j) ) 50 continue * reflection boundary condition. call reflec(pf,nx,ny) end subroutine projec(uf,vf,us,vs,px,py,nx,ny,mxy,piter) c******************************************************************* c c correct the velocity field c c (uf,vf) <- generated velocity field c (us,vs) <- velocity correction before project c grad cor <- gradient component of (us,vs) c c local arrays: c cor,res,div: 1-d arrays used in the multigrid package c to store the solution, residue and the rhs. c (adx,ady): grad of cor. c c******************************************************************* integer nx,ny,i,j,k,kk,nxm,nym,mxy,mxym include 'gridsz' logical piter double precision uf(nx,ny),vf(nx,ny), + us(nx,ny),vs(nx,ny), + px(nx,ny),py(nx,ny) double precision cor(mxym),res(mxym),div(mxym),ro(mxym) double precision adx(nxm,nym),ady(nxm,nym) double precision hhinv,dt,roa common/timest/dt common/hfract/hhinv common/density/roa(0:nxm+1,0:nym+1) data cor/mxym*0.D0/ save cor k(i,j) = (nx+3)*(j+1)+i+2 do 2 i=1,mxym div(i) = 0.D0 res(i) = 0.D0 2 continue kk = nx+5 do 3 j=0,ny+1 do 4 i=0,nx+1 ro(kk) = 1.D0/roa(i,j) kk = kk+1 4 continue kk = kk+1 3 continue c*********************************************************** c obtain the field to be projected do 5 j=1,ny do 5 i=1,nx us(i,j) = (us(i,j)-uf(i,j))/dt vs(i,j) = (vs(i,j)-vf(i,j))/dt 5 continue c call formdv(us,vs,div,nx,ny) c call multi2(div,cor,ro,res,nx,ny,mxy) c c generate the gradient of cor do 10 j=1,ny do 10 i=1,nx adx(i,j) = (cor(k(i,j)) - cor(k(i-1,j)) + + cor(k(i,j-1)) - cor(k(i-1,j-1)) )*hhinv ady(i,j) = (cor(k(i,j)) - cor(k(i,j-1)) + + cor(k(i-1,j)) - cor(k(i-1,j-1)) )*hhinv 10 continue c*************************************************** c update the pressure do 30 j=1,ny do 30 i=1,nx px(i,j) = px(i,j) + adx(i,j) py(i,j) = py(i,j) + ady(i,j) 30 continue c c update the velocity if (.not.piter) then do 40 j=1,ny do 40 i=1,nx uf(i,j) = uf(i,j) + (us(i,j) - adx(i,j)/roa(i,j))*dt vf(i,j) = vf(i,j) + (vs(i,j) - ady(i,j)/roa(i,j))*dt 40 continue end if end subroutine formdv(u,v,div,nx,ny) c******************************************************************** c c (u,v) <- velocity field c div -> h^2 div (u,v) c assuming homogeneous bc. c c******************************************************************** integer i,j,nx,ny double precision u(nx,ny),v(nx,ny), + div(-1:nx+1,-1:ny+1) double precision h,hinv,h2,sdiv,divm,dt common/timest/dt common/meshsz/h,hinv h2 = 0.5D0*h do 10 j=1,ny-1 do 10 i=1,nx-1 div(i,j) = h2*( + u(i+1,j+1) - u(i,j+1) + + u(i+1,j ) - u(i,j ) + + v(i+1,j+1) - v(i+1,j) + + v(i, j+1) - v(i, j) ) 10 continue do 20 j=1,ny-1 div(0,j) = h2*( + u(1,j+1) + u(1,j) + + v(1,j+1) - v(1,j) ) div(nx,j) = h2*( + -u(nx,j+1) - u(nx,j) + + v(nx,j+1) - v(nx,j) ) 20 continue do 30 i=1,nx-1 div(i,0) = h2*( + u(i+1,1) - u(i,1) + + v(i+1,1) + v(i,1) ) div(i,ny) = h2*( + u(i+1,ny) - u(i,ny) - + v(i+1,ny) - v(i,ny)) 30 continue div(0,0) = h2*(u(1,1) + v(1,1)) div(0,ny) = h2*(u(1,ny) - v(1,ny)) div(nx,0) = -h2*(u(nx,1) - v(nx,1)) div(nx,ny) = -h2*(u(nx,ny) + v(nx,ny)) c sdiv = 0.D0 divm = 0.D0 do 110 j=1,ny do 110 i=1,nx divm = dmax1(divm,dabs(div(i,j))) sdiv = sdiv + div(i,j)*div(i,j) 110 continue sdiv = dsqrt(sdiv)*hinv divm = divm*hinv*hinv print *,'l2 of div',sdiv,' divm=',divm end subroutine macpro(ux,vx,uy,vy,nx,ny,mxy) ********************************************************** * correcting the edge velocities by subtracting the * pressure gradient from the mac projection. * * ---(uy,vy)---- * | * | * (i,j) (ux,vx) * | * | * -------------- *********************************************************** integer i,j,nx,ny,k,nxm,nym,mxym,mxy include 'gridsz' double precision ux(0:nx,0:ny),vx(0:nx,0:ny), + uy(0:nx,0:ny),vy(0:nx,0:ny) double precision adx(0:nxm,nym),ady(nxm,0:nym), + cor(mxym),res(mxym),div(mxym),ro(mxym) double precision h,hinv,roa,rr,rt common/meshsz/h,hinv common/density/roa(0:nxm+1,0:nym+1) data cor/mxym*0.D0/ save cor k(i,j) = (nx+2)*j+i+1 do 2 i=1,mxym res(i) = 0.D0 div(i) = 0.D0 2 continue do 3 j=0,ny+1 do 3 i=0,nx+1 ro(k(i,j)) = 1.D0/roa(i,j) 3 continue c obtain the divergence for the edge velocity field (MAC projection), c multiplied by h^2 call fmdved(ux,vx,uy,vy,div,nx,ny) call multi1(div,cor,ro,res,nx,ny,mxy) do 10 j=1,ny do 10 i=1,nx-1 rr = 0.5D0*(ro(k(i,j))+ro(k(i+1,j))) adx(i,j) = hinv*(cor(k(i+1,j))-cor(k(i,j)))*rr 10 continue do 12 j=1,ny adx(0,j) = 0.D0 adx(nx,j) = 0.D0 12 continue do 20 j=1,ny-1 do 20 i=1,nx rt = 0.5D0*(ro(k(i,j))+ro(k(i,j+1))) ady(i,j) = hinv*(cor(k(i,j+1))-cor(k(i,j)))*rt 20 continue do 22 i=1,nx ady(i,0) = 0.D0 ady(i,ny) = 0.D0 22 continue do 30 j=1,ny do 30 i=1,nx-1 ux(i,j) = ux(i,j) - adx(i,j) 30 continue do 40 j=1,ny-1 do 40 i=1,nx vy(i,j) = vy(i,j) - ady(i,j) 40 continue do 50 j=1,ny-1 do 50 i=1,nx uy(i,j) = uy(i,j) - 0.25D0*( + adx(i,j) + adx(i-1,j) + + adx(i,j+1) + adx(i-1,j+1) ) 50 continue do 60 j=1,ny do 60 i=1,nx-1 vx(i,j) = vx(i,j) - 0.25D0*( + ady(i,j) + ady(i+1,j) + + ady(i,j-1) + ady(i+1,j-1) ) 60 continue end subroutine fmdved(ux,vx,uy,vy,div,nx,ny) integer i,j,nx,ny logical prnorm double precision ux(0:nx,0:ny),vx(0:nx,0:ny), + uy(0:nx,0:ny),vy(0:nx,0:ny), + div(0:nx+1,0:ny+1) double precision divm,sdiv,h,hinv common/meshsz/h,hinv parameter (prnorm=.false.) do 5 j=1,ny do 5 i=1,nx div(i,j) = h*( ux(i,j) - ux(i-1,j) + + vy(i,j) - vy(i,j-1) ) 5 continue C if (prnorm .eq. .false.) return if (.not. prnorm) return divm = 0.D0 sdiv = 0.D0 do 10 j=1,ny do 10 i=1,nx divm = dmax1(divm,dabs(div(i,j))) sdiv = sdiv + div(i,j)**2 10 continue sdiv = dsqrt(sdiv)*hinv divm = divm*hinv*hinv print *,'mac l2=',sdiv,' l_inf=',divm end subroutine reinit(phi,nx,ny) ******************************************************************* * reinitialize the level set function so it is the distance * function from the zero level set. * Given phi on the grid, calculate the time it takes for a * propagating front with propagating speed 1, which starts as * the zero level set of phi, to arrive at the grid point (i,j). * The time for this particular grid point is stored in tarv(i,j). * Two sides of the zero level set are dealt with separately in * two sweeps (k=1,2 in loop 2). * * phi: level set function; * fo,fn: local arrays used to contain the propagating * information of the "front", fo (initial) and * fn (final) states. ******************************************************************* integer nx,ny,nxm,nym,mxym,i,j,k,ndone,kdel,kt,kk include 'gridsz' parameter(kdel=10) logical done(0:nxm+1,0:nym+1) double precision t,dt,del,h,hinv,sign,dth double precision phi(0:nx+1,0:ny+1) double precision fo(0:nxm+1,0:nym+1),fn(0:nxm+1,0:nym+1), + tarv(0:nxm+1,0:nym+1) common/meshsz/h,hinv del = dble(kdel)*h dt = 0.5D0*h dth = 0.5D0 kt = 2*kdel c do 2 kk=1,2 sign = (-1.D0)**(kk+1) ndone = 0 t = 0.D0 do 8 j=0,ny+1 do 8 i=0,nx+1 fo(i,j) = sign*phi(i,j) 8 continue c mark all the points on the nonnegative side: do 10 j=0,ny+1 do 10 i=0,nx+1 done(i,j) = fo(i,j) .ge. 0.D0 if(done(i,j)) ndone = ndone + 1 10 continue * do 5 k=1,kt * do 50 j=1,ny do 50 i=1,nx fn(i,j) = fo(i,j) + dth*dsqrt( + dmin1(fo(i,j)-fo(i-1,j),0.D0)**2 + + dmax1(fo(i+1,j)-fo(i,j),0.D0)**2 + + dmin1(fo(i,j)-fo(i,j-1),0.D0)**2 + + dmax1(fo(i,j+1)-fo(i,j),0.D0)**2 ) 50 continue do 60 j=1,ny fn(0,j) = fn(1,j) -h fn(nx+1,j) = fn(nx,j) - h 60 continue do 65 i=1,nx fn(i,0) = fn(i,1) - h fn(i,ny+1) = fn(i,ny) - h 65 continue fn(0,0) = fn(1,1) - h fn(nx+1,0) = fn(nx,1) - h fn(0,ny+1) = fn(1,ny) - h fn(nx+1,ny+1) = fn(nx,ny) - h c t = t + dt c do 20 j=0,ny+1 do 20 i=0,nx+1 if( .not. done(i,j) .and. fn(i,j) .ge. 0.D0 ) then tarv(i,j) = t - dt*fn(i,j)/(fn(i,j)-fo(i,j)) ndone = ndone + 1 done(i,j) = .true. endif 20 continue c if(ndone .eq. (nx+2)*(ny+2)) goto 2 c do 25 j=0,ny+1 do 25 i=0,nx+1 fo(i,j) = fn(i,j) 25 continue c 5 continue * do 30 j=0,ny+1 do 30 i=0,nx+1 if( .not. done(i,j) ) tarv(i,j) = del 30 continue 2 continue do 70 j=0,ny+1 do 70 i=0,nx+1 phi(i,j) = dsign(tarv(i,j),phi(i,j)) 70 continue end subroutine vortic(u,v,vor,nx,ny) c******************************************************************** c c obtain the vorticity of the velocity field. c c (u,v) -> given velocity field c vor -> vorticity field. c c******************************************************************** integer i,j,nx,ny double precision u(nx,ny),v(nx,ny),vor(nx,ny) double precision hhinv,thinv common/hfract/hhinv thinv = 2.D0/3.D0*hhinv do 10 j=2,ny-1 do 10 i=2,nx-1 vor(i,j) = hhinv*( v(i+1,j) - v(i-1,j) - + u(i,j+1) + u(i,j-1)) 10 continue do 20 j=2,ny-1 vor(1,j) = thinv*( v(2,j) + 3.D0*v(1,j) ) + + hhinv*( u(1,j-1) - u(1,j+1) ) vor(nx,j) = -thinv*( v(nx-1,j) + 3.D0*v(nx,j) ) + + hhinv*( u(nx,j-1) - u(nx,j+1) ) 20 continue do 30 i=2,nx-1 vor(i,1) = hhinv*( v(i+1,1) - v(i-1,1) ) - + thinv*( u(i,2) + 3.D0*u(i,1) ) vor(i,ny) = hhinv*( v(i+1,ny) - v(i-1,ny) ) + + thinv*( u(i,ny-1) + 3.D0*u(i,ny) ) 30 continue vor(1,1) = thinv*( v(2,1) + 3.D0*v(1,1) - + u(1,2) - 3.D0*u(1,1) ) vor(1,ny) = thinv*( v(2,ny) + 3.D0*v(1,ny) + + u(1,ny-1) + 3.D0*u(1,ny) ) vor(nx,1) = thinv*( -v(nx-1,1) - 3.D0*v(nx,1) - + u(nx,2) - 3.D0*u(nx,1) ) vor(nx,ny) = thinv*( -v(nx-1,ny) - 3.D0*v(nx,ny) + + u(nx,ny-1) + 3.D0*u(nx,ny) ) c write(24,*) vor end subroutine chtmst(u,v,ek,dtin,cfl,nx,ny) c***************************************************************** c c impose the cfl condition on the time step and calculate the c total kinetic energy. c c dtin -> maximal time step c cfl -> cfl number c dt <- time step c ek <- total kinetic energy c c***************************************************************** integer i,j,nx,ny double precision u(nx,ny),v(nx,ny), + h,hinv,ek,dtin,cfl,uvmax,dt common/meshsz/h,hinv common/timest/dt uvmax = 0.D0 ek = 0.D0 do 10 j=1,ny do 10 i=1,nx uvmax = dmax1(dabs(u(i,j)),dabs(v(i,j)),uvmax) ek = ek + u(i,j)*u(i,j) + v(i,j)*v(i,j) 10 continue * ek = 0.5D0*ek*h*h * total kinetic energy if(uvmax .eq. 0.D0) then dt = dtin else dt = dmin1(dtin,cfl*h/uvmax) endif end subroutine output(pf,u,v,nx,ny,kunit) integer i,j,nx,ny,kunit double precision pf(0:nx+1,0:ny+1),u(nx,ny),v(nx,ny) c write(21,101) ((u(i,j),i=1,nx),j=1,ny) c write(21,101) ((v(i,j),i=1,nx),j=1,ny) print *,'kunit=',kunit write(kunit,*) pf 101 format(e16.10) 102 format(2(e16.10,2x)) 103 format(/) end subroutine upwind(ul,ur,vl,vr,u,v) c **************************************************************** c c Godunov upwinding c c **************************************************************** double precision ul,ur,vl,vr,u,v if (ul .ge. 0.D+00 .and. ul+ur .ge. 0.D+00) then u = ul elseif ( ul .le. 0.D+00 .and. ur .ge. 0.D+00) then u = 0.D+00 else u = ur endif if (u) 10,20,30 10 v = vr return 20 v = 0.5D+00*(vl+vr) return 30 v = vl return end subroutine upwsec(u,ul,ur,v,vl,vr,du,dul,dur,dv,dvl,dvr, + duup,dvup) c **************************************************************** c c second order upwinding of the derivative c c **************************************************************** double precision u,ul,ur,v,vl,vr,du,dul,dur,dv,dvl,dvr, + duup,dvup,dth,duc common/timspc/dth if (u .ge. 0.D0) then duc = 0.5D0*( 1.D0 - u*dth ) duup = u - ul + duc*( du - dul ) dvup = v - vl + duc*( dv - dvl ) else duc = -0.5D0*( 1.D0 + u*dth ) duup = ur - u + duc*( dur - du ) dvup = vr - v + duc*( dvr - dv ) endif end function dm(dr,dl) * ******************************************************** * * limiter for the second-order hj upwind scheme * * ******************************************************** double precision dm,dr,dl if( dabs(dr) .le. dabs(dl) ) then dm = dr else dm = dl endif end function slimit(ul,ur) c *********************************************************** c c limiter definition in the interior c c ul,ur -> values on the left and right sides c slimit <- resolved value c c *********************************************************** double precision ul,ur,dut,slimit if ( ul*ur .le. 0.D+00 ) then slimit = 0.D+00 else dut = 0.5D0*( ul + ur ) slimit = dsign( 1.D0, dut )*dmin1( + dabs(dut), 2.D0*dabs(ul), 2.D0*dabs(ur) ) end if c slimit = 0.5D0*(ul+ur) end function blimit(ub,u1,u2) c *********************************************************** c c limiter definition on the boundaries c c ub,u1,u2 -> values spaced by h/2,h and h c blimit <- resolved value c c *********************************************************** double precision ub,u1,u2,blimit,umax,umin,ubt ubt = 0.25D0*( 3.D0*u1 - u2 + 2.D0*ub ) umax = dmax1( u1, ub ) umin = dmin1( u1, ub ) if( ubt .gt. umax ) then blimit = 2.D0*( u1 - umax ) else if( ubt .lt. umin ) then blimit = 2.D0*( u1 - umin ) else blimit = 0.5D0*( u2 + u1 - 2.D0*ub ) end if end subroutine inifrt(pf,rin,nx,ny) ***************************************************** c initialize the level set function ***************************************************** integer i,j,nx,ny double precision pf(0:nx+1,0:ny+1) double precision h,hinv,x,y,r,rin common/meshsz/h,hinv do 10 j=0,ny+1 y = (dfloat(j)-0.5D0)*h do 10 i=0,nx+1 x = (dfloat(i)-0.5D0)*h r = dsqrt( (x-0.5D0)**2 + (y-0.5D0)**2 ) pf(i,j) = rin - r 10 continue end subroutine reflec(pf,nx,ny) * *********************************************** * impose the reflection boundary condition for pf * *********************************************** integer i,j,nx,ny double precision pf(0:nx+1,0:ny+1) do 10 i=1,nx pf(i,0) = pf(i,1) pf(i,ny+1) = pf(i,ny) 10 continue do 20 j=1,ny pf(0,j) = pf(1,j) pf(nx+1,j) = pf(nx,j) 20 continue end subroutine multi1(f,u,ro,res,nx,ny,mxy) ***************************************************************** * * MULTIGRID PACKAGE TO SOLVE THE 5-POINT VARIABLE COEFFECIENT * LAPLACIAN: * * div ( ro grad u ) = f * * USED IN MAC PROJECTION, where res is the residue. Simple V * cycles are used. * * parameters: * kit: number of iterations at each level, * nv: max number of V cycles. ***************************************************************** integer mxy,nx,ny,nnx,nny,nlevel,j,lnx,lny,nv,mnx,mny, + kf,kc,kit,kv,nxy2 integer kp(10) double precision f(mxy),u(mxy),ro(mxy),res(mxy) double precision el2,h,hinv,hh,eps common/meshsz/h,hinv parameter (kit=2,nv=8) parameter (eps=1.D-10) call dumpnu(f,nx,ny) nnx = nx nny = ny nlevel = 1 kp(1) = 1 * organize the space storage and obtain the pointers do 20 j=2,10 if (nny .lt. 2) goto 15 nlevel = nlevel + 1 kp(j) = kp(j-1) + (nnx+2)*(nny+2) call resde1(nnx,nny,nnx/2,nny/2,ro(kp(j-1)),ro(kp(j))) nnx = nnx/2 nny = nny/2 20 continue 15 nnx = nx nny = ny do 5 kv=1,nv hh = h c descending step: do 25 j=1,nlevel-1 kf = kp(j) kc = kp(j+1) nxy2 = (nnx+2)*(nny+2) call gsneu1(nnx,nny,nxy2,u(kf),f(kf),ro(kf),kit) call resid1(nnx,nny,nxy2,hh,f(kf),u(kf),ro(kf),res(kf),el2) c print *,'j=',j,' el2=',el2 lnx = nnx/2 lny = nny/2 call restr1(nnx,nny,lnx,lny,res(kf),f(kc)) nnx = lny nny = lny hh = 2.D0*hh 25 continue call dirso1(nnx,nny,u(kc),f(kc)) do 26 j=nlevel,2,-1 kc = kp(j) kf = kp(j-1) mnx = 2*nnx mny = 2*nny call intpl1(nnx,nny,mnx,mny,u(kc),u(kf)) nxy2 = (nnx+2)*(nny+2) call setzro(u(kc),nxy2) nnx = mnx nny = mny nxy2 = (nnx+2)*(nny+2) call gsneu1(nnx,nny,nxy2,u(kf),f(kf),ro(kf),kit) 26 continue call resid1(nx,ny,nxy2,h,f,u,ro,res,el2) c print *,'cycle ',kv,' el2= ',el2 if (el2 .le. eps) goto 6 call dumpnu(u,nx,ny) 5 continue 6 print *,'projection el2',el2 101 format(2(f12.6,2x)) end subroutine gsneu1(nx,ny,nxy2,u,f,ro,kit) **************************************************************** * red-black G-S iterations, with Neumann boundary conditions. **************************************************************** integer nx,ny,i,j,k,kv,kit,nx2,nxy2 double precision f(nxy2),u(nxy2),ro(nxy2) double precision rr,rl,rb,rt,rs nx2 = nx+2 do 5 kv=1,kit call presb1(u,nx,ny,nxy2) do 30 j=1,ny-1,2 do 30 i=1,nx-1,2 k = nx2*j+i+1 rr = 0.5D0*(ro(k) + ro(k+1)) rl = 0.5D0*(ro(k) + ro(k-1)) rb = 0.5D0*(ro(k) + ro(k-nx2)) rt = 0.5D0*(ro(k) + ro(k+nx2)) rs = rr + rl + rb + rt u(k) = (rr*u(k+1) + rl*u(k-1) + + rt*u(k+nx2) + rb*u(k-nx2) - + f(k))/rs 30 continue do 40 j=2,ny,2 do 40 i=2,nx,2 k = nx2*j+i+1 rr = 0.5D0*(ro(k) + ro(k+1)) rl = 0.5D0*(ro(k) + ro(k-1)) rb = 0.5D0*(ro(k) + ro(k-nx2)) rt = 0.5D0*(ro(k) + ro(k+nx2)) rs = rr + rl + rb + rt u(k) = (rr*u(k+1) + rl*u(k-1) + + rt*u(k+nx2) + rb*u(k-nx2) - + f(k))/rs 40 continue call presb1(u,nx,ny,nxy2) do 50 j=1,ny-1,2 do 50 i=2,nx,2 k = nx2*j+i+1 rr = 0.5D0*(ro(k) + ro(k+1)) rl = 0.5D0*(ro(k) + ro(k-1)) rb = 0.5D0*(ro(k) + ro(k-nx2)) rt = 0.5D0*(ro(k) + ro(k+nx2)) rs = rr + rl + rb + rt u(k) = (rr*u(k+1) + rl*u(k-1) + + rt*u(k+nx2) + rb*u(k-nx2) - + f(k))/rs 50 continue do 60 j=2,ny,2 do 60 i=1,nx-1,2 k = nx2*j+i+1 rr = 0.5D0*(ro(k) + ro(k+1)) rl = 0.5D0*(ro(k) + ro(k-1)) rb = 0.5D0*(ro(k) + ro(k-nx2)) rt = 0.5D0*(ro(k) + ro(k+nx2)) rs = rr + rl + rb + rt u(k) = (rr*u(k+1) + rl*u(k-1) + + rt*u(k+nx2) + rb*u(k-nx2) - + f(k))/rs 60 continue 5 continue end subroutine resid1(nx,ny,nxy2,h,f,u,ro,res,el2) ****************************************************************** * compute the residue ******************************************************* integer nx,ny,i,j,k,nx2,nxy2 double precision f(nxy2),res(nxy2),u(nxy2),ro(nxy2) double precision el2,h,rr,rl,rb,rt,rs nx2 = nx+2 call presb1(u,nx,ny,nxy2) do 10 j=1,ny do 12 i=1,nx k = nx2*j+i+1 rr = 0.5D0*(ro(k) + ro(k+1)) rl = 0.5D0*(ro(k) + ro(k-1)) rb = 0.5D0*(ro(k) + ro(k-nx2)) rt = 0.5D0*(ro(k) + ro(k+nx2)) rs = rr + rl + rb + rt res(k) = f(k) + rs*u(k) - + rr*u(k+1) - rl*u(k-1) - + rt*u(k+nx2) - rb*u(k-nx2) 12 continue 10 continue el2 = 0.D0 do 60 k=1,nx*ny el2 = el2 + res(k)*res(k) 60 continue el2 = dsqrt(el2)/h end subroutine intpl1(nx,ny,mx,my,uc,uf) ******************************************************** * interpolate from the coarse to fine grid * mx=2*nx, my=2*ny ******************************************************** integer nx,ny,mx,my,i,j,kc,kf,nx2,mx2 double precision uc((nx+2)*(ny+2)),uf((mx+2)*(my+2)) double precision uadd nx2 = nx+2 mx2 = mx+2 kc = nx+4 kf = mx+4 do 10 j=1,ny do 11 i=1,nx uadd = uc(kc) uf(kf) = uf(kf) + uadd uf(kf+1) = uf(kf+1) + uadd uf(kf+mx2) = uf(kf+mx2) + uadd uf(kf+mx2+1) = uf(kf+mx2+1) + uadd kc = kc+1 kf = kf+2 11 continue kc = kc+2 kf = kf+mx2+2 10 continue end subroutine restr1(nx,ny,lx,ly,rf,rc) *********************************************************** * restrict from the fine to coarse grid * lx=nx/2, ly=ny/2. *********************************************************** integer nx,ny,lx,ly,nx2,lx2,i,j,kc,kf double precision rf((nx+2)*(ny+2)),rc((lx+2)*(ly+2)) nx2 = nx+2 lx2 = lx+2 kc = lx+4 kf = nx+4 do 10 j=1,ly do 11 i=1,lx rc(kc) = rf(kf) + rf(kf+1) + + rf(kf+nx2) + rf(kf+nx2+1) kc = kc+1 kf = kf+2 11 continue kc = kc+2 kf = kf+nx2+2 10 continue end subroutine resde1(nx,ny,lx,ly,df,dc) ************************************************************ * averaging the density df in 4 neighboring cells in the fine * grid to obtain the density dc in the coarse grid. * lx = nx/2, ly=ny/2 ************************************************************ integer nx,ny,lx,ly,i,j,kc,kf,kc1,kc2,nx2,lx2 double precision df((nx+2)*(ny+2)),dc((lx+2)*(ly+2)) nx2 = nx+2 lx2 = lx+2 kc = lx+4 kf = nx+4 do 10 j=1,ly do 11 i=1,lx dc(kc) = 0.25D0*(df(kf) + df(kf+1) + + df(kf+nx2) + df(kf+nx2+1) ) kc = kc+1 kf = kf+2 11 continue kc = kc+2 kf = kf+nx2+2 10 continue kc1 = 2 kc2 = lx2*(ly+1)+2 do 20 i=1,lx dc(kc1) = dc(kc1+lx2) dc(kc2) = dc(kc2-lx2) kc1 = kc1+1 kc2 = kc2+1 20 continue kc1 = 1 kc2 = lx2 do 30 j=0,ly+1 dc(kc1) = dc(kc1+1) dc(kc2) = dc(kc2-1) kc1 = kc1+lx2 kc2 = kc2+lx2 30 continue end subroutine dirso1(nx,ny,u,f) ************************************************************ * solving the coarse grid exactly by inverting the matrix. ************************************************************ integer nx,ny double precision u((nx+2)*(ny+2)),f(nx*ny) u(nx+4) = 0.D0 end subroutine presb1(u,nx,ny,nxy2) ************************************************************* * prescirbe the Neumann boundary condition ************************************************************* integer k1,k2,i,j,nx,ny,nxy2,nx2 double precision u(nxy2) nx2 = nx+2 k1 = 1 k2 = nx2*(ny+1)+1 do 10 i=1,nx k1 = k1+1 u(k1) = u(k1+nx2) k2 = k2+1 u(k2) = u(k2-nx2) 10 continue k1 = 1 k2 = nx2 do 20 j=1,ny k1 = k1+nx2 u(k1) = u(k1+1) k2 = k2+nx2 u(k2) = u(k2-1) 20 continue end subroutine dumpnu(f,nx,ny) ************************************************************ c force f to be mean zero, by subtract its average. ************************************************************ integer k,i,j,nx,ny double precision f(0:nx+1,0:ny+1),fsum fsum = 0.D0 k = 0 do 10 j=1,ny do 10 i=1,nx k = k+1 fsum = fsum + (f(i,j)-fsum)/dble(k) 10 continue do 20 j=1,ny do 20 i=1,nx f(i,j) = f(i,j) - fsum 20 continue end subroutine multi2(f,u,ro,res,nx,ny,mxy) ****************************************************************** * * MULTIGRID PACKAGE FOR SOLVING 5-POINT LAPLACIAN * USED IN THE APPROXIMATE PROJECTION STEP. * ****************************************************************** integer mxy,nx,ny,nnx,nny,nlevel,j,lnx,lny,nv,mnx,mny, + kf,kc,kit,kv,nxy3 integer kp(10) logical prnorm double precision f(mxy),u(mxy),ro(mxy),res(mxy) double precision el2,h,hinv,hh,eps common/meshsz/h,hinv parameter (kit=2,nv=8) parameter (eps=1.D-10) parameter (prnorm=.true.) call dumpf2(f,nx,ny) nnx = nx nny = ny nlevel = 1 kp(1) = 1 * organize the space storage and obtain the pointers do 20 j=2,10 if (nny .lt. 2) goto 15 nlevel = nlevel + 1 kp(j) = kp(j-1) + (nnx+3)*(nny+3) call resden(nnx,nny,nnx/2,nny/2,ro(kp(j-1)),ro(kp(j))) nnx = nnx/2 nny = nny/2 20 continue 15 nnx = nx nny = ny do 5 kv=1,nv hh = h * descending steps: do 25 j=1,nlevel-1 kf = kp(j) kc = kp(j+1) nxy3 = (nnx+3)*(nny+3) call gsneu2(nnx,nny,nxy3,u(kf),f(kf),ro(kf),kit) call resid2(nnx,nny,nxy3,hh,f(kf),u(kf),ro(kf),res(kf),el2) c print *,'j=',j,' el2=',el2 lnx = nnx/2 lny = nny/2 call restr2(nnx,nny,lnx,lny,res(kf),f(kc)) nnx = lny nny = lny hh = 2.D0*hh 25 continue * solve the coarse grid equations exactly by inverting the * matrix nxy3 = (nnx+3)*(nny+3) call dirso2(nnx,nny,nxy3,u(kc),f(kc),ro(kc)) * ascending step: do 26 j=nlevel,2,-1 kc = kp(j) kf = kp(j-1) mnx = 2*nnx mny = 2*nny call intpl2(nnx,nny,mnx,mny,u(kc),u(kf)) nxy3 = (nnx+3)*(nny+3) call setzro(u(kc),nxy3) nnx = mnx nny = mny nxy3 = (nnx+3)*(nny+3) call gsneu2(nnx,nny,nxy3,u(kf),f(kf),ro(kf),kit) 26 continue call resid2(nx,ny,nxy3,h,f,u,ro,res,el2) c print *,'cycle ',kv,' el2= ',el2 if (el2 .le. eps) goto 6 call dumpu2(u,nx,ny) 5 continue C6 if(prnorm .eq. .true.) print *,'projection el2',el2 6 if(prnorm) print *,'projection el2',el2 101 format(2(f12.6,2x)) end subroutine gsneu2(nx,ny,nxy3,u,f,ro,kit) ***************************************************************** * red-black G-S iterations, with Neumann boundary conditions. ***************************************************************** integer nx,ny,nxy3,i,j,k,kp,kv,kit,nx3 double precision f(nxy3),u(nxy3),ro(nxy3) double precision rr,rl,rb,rt,rs kp(i,j) = (nx+3)*(j+1)+i+2 nx3 = nx+3 do 5 kv=1,kit * prescribe the Neumann boundary condition: call presbd(u,ro,nx,ny,nxy3) do 30 j=0,ny,2 do 30 i=0,nx,2 k = kp(i,j) rr = 0.5D0*(ro(k+1) + ro(k+nx3+1)) rl = 0.5D0*(ro(k) + ro(k+nx3)) rb = 0.5D0*(ro(k) + ro(k+1)) rt = 0.5D0*(ro(k+nx3) + ro(k+nx3+1)) rs = ro(k) + ro(k+1) + ro(k+nx3) + ro(k+nx3+1) u(k) = (rr*u(k+1) + rl*u(k-1) + + rt*u(k+nx3) + rb*u(k-nx3) - + f(k) )/rs 30 continue do 40 j=1,ny-1,2 do 40 i=1,nx-1,2 k = kp(i,j) rr = 0.5D0*(ro(k+1) + ro(k+nx3+1)) rl = 0.5D0*(ro(k) + ro(k+nx3)) rb = 0.5D0*(ro(k) + ro(k+1)) rt = 0.5D0*(ro(k+nx3) + ro(k+nx3+1)) rs = ro(k) + ro(k+1) + ro(k+nx3) + ro(k+nx3+1) u(k) = (rr*u(k+1) + rl*u(k-1) + + rt*u(k+nx3) + rb*u(k-nx3) - + f(k) )/rs 40 continue call presbd(u,ro,nx,ny,nxy3) do 50 j=0,ny,2 do 50 i=1,nx-1,2 k = kp(i,j) rr = 0.5D0*(ro(k+1) + ro(k+nx3+1)) rl = 0.5D0*(ro(k) + ro(k+nx3)) rb = 0.5D0*(ro(k) + ro(k+1)) rt = 0.5D0*(ro(k+nx3) + ro(k+nx3+1)) rs = ro(k) + ro(k+1) + ro(k+nx3) + ro(k+nx3+1) u(k) = (rr*u(k+1) + rl*u(k-1) + + rt*u(k+nx3) + rb*u(k-nx3) - + f(k) )/rs 50 continue do 60 j=1,ny-1,2 do 60 i=0,nx,2 k = kp(i,j) rr = 0.5D0*(ro(k+1) + ro(k+nx3+1)) rl = 0.5D0*(ro(k) + ro(k+nx3)) rb = 0.5D0*(ro(k) + ro(k+1)) rt = 0.5D0*(ro(k+nx3) + ro(k+nx3+1)) rs = ro(k) + ro(k+1) + ro(k+nx3) + ro(k+nx3+1) u(k) = (rr*u(k+1) + rl*u(k-1) + + rt*u(k+nx3) + rb*u(k-nx3) - + f(k) )/rs 60 continue 5 continue end subroutine resid2(nx,ny,nxy3,h,f,u,ro,res,el2) ******************************************************* * compute the residue ******************************************************* integer nx,ny,i,j,k,nx3,nxy3 double precision f(nxy3),res(nxy3),u(nxy3),ro(nxy3) double precision el2,h,rr,rl,rb,rt,rs nx3 = nx+3 call presbd(u,ro,nx,ny,nxy3) el2 = 0.D0 do 10 j=0,ny do 12 i=0,nx k = nx3*(j+1)+i+2 rr = 0.5D0*(ro(k+1) + ro(k+nx3+1)) rl = 0.5D0*(ro(k) + ro(k+nx3)) rb = 0.5D0*(ro(k) + ro(k+1)) rt = 0.5D0*(ro(k+nx3) + ro(k+nx3+1)) rs = ro(k) + ro(k+1) + ro(k+nx3) + ro(k+nx3+1) res(k) = f(k) + rs*u(k) - + rr*u(k+1) - rl*u(k-1) - + rt*u(k+nx3) - rb*u(k-nx3) el2 = el2 + res(k)*res(k) 12 continue 10 continue c el2 = sqrt( (res(i,j)/h^2)^2 *h^2 ) el2 = dsqrt(el2)/h end subroutine setzro(u,nxy) integer k,nxy double precision u(nxy) do 10 k=1,nxy u(k) = 0.D0 10 continue end subroutine intpl2(nx,ny,mx,my,uc,uf) ******************************************************** * interpolate from the coarse to fine grid * mx=2*nx, my=2*ny ******************************************************** integer nx,ny,mx,my,i,j,kc,kf,mx2,mx3,mx4 double precision uc(*),uf(*) double precision uadd mx3 = mx+3 mx2 = mx+2 mx4 = mx+4 kc=nx+5 kf=mx+5 do 10 j=0,ny do 11 i=0,nx uadd = uc(kc) uf(kf) = uf(kf) + uadd uf(kf+1) = uf(kf+1) + 0.5D0*uadd uf(kf-1) = uf(kf-1) + 0.5D0*uadd uf(kf+mx3) = uf(kf+mx3) + 0.5D0*uadd uf(kf+mx2) = uf(kf+mx2) + 0.25D0*uadd uf(kf+mx4) = uf(kf+mx4) + 0.25D0*uadd uf(kf-mx3) = uf(kf-mx3) + 0.5D0*uadd uf(kf-mx2) = uf(kf-mx2) + 0.25D0*uadd uf(kf-mx4) = uf(kf-mx4) + 0.25D0*uadd kc = kc + 1 kf = kf + 2 11 continue kc = kc + 2 kf = kf + mx3 + 1 10 continue end subroutine restr2(nx,ny,lx,ly,rf,rc) *********************************************************** * restrict from the fine to coarse grid * lx=nx/2, ly=ny/2. *********************************************************** integer nx,ny,lx,ly,i,j,kc,kf,nx3 double precision rf(*),rc(*) nx3 = nx+3 kc = lx+5 kf = nx+5 do 10 j=0,ly do 11 i=0,lx rc(kc) = rf(kf) + 0.5D0*( + rf(kf+1) + rf(kf-1) + + rf(kf+nx3) + rf(kf-nx3) ) + 0.25D0*( + rf(kf+nx3+1) + rf(kf+nx3-1) + + rf(kf-nx3+1) + rf(kf-nx3-1) ) kc = kc + 1 kf = kf + 2 11 continue kc = kc + 2 kf = kf + nx3 + 1 10 continue end subroutine resden(nx,ny,lx,ly,df,dc) ************************************************************ * averaging the density df in 4 neighboring cells in the fine * grid to obtain the density dc in the coarse grid. * lx = nx/2, ly=ny/2 ************************************************************ integer nx,ny,lx,ly,i,j,kc,kf,kc1,kc2,nx3,lx3 double precision df(*),dc(*) nx3 = nx+3 lx3 = lx+3 kc = 2*lx3+3 kf = 2*nx3+3 do 10 j=1,ly do 11 i=1,lx dc(kc) = 0.25D0*( df(kf) + df(kf+1) + + df(kf+nx3) + df(kf+nx3+1) ) kc = kc + 1 kf = kf + 2 11 continue kc = kc + 3 kf = kf + nx3 + 3 10 continue kc1 = lx3+3 kc2 = lx3*(ly+2)+3 do 20 i=1,lx dc(kc1) = dc(kc1+lx3) dc(kc2) = dc(kc2-lx3) kc1 = kc1+1 kc2 = kc2+1 20 continue kc1 = lx3+2 kc2 = 2*lx3 do 30 j=0,ly+1 dc(kc1) = dc(kc1+1) dc(kc2) = dc(kc2-1) kc1 = kc1+lx3 kc2 = kc2+lx3 30 continue end subroutine dirso2(nx,ny,nxy3,u,f,ro) ************************************************************ * solving the coarse grid exactly by inverting the matrix. ************************************************************ integer nx,ny,nxy3,k double precision u(nxy3),f(nxy3),ro(nxy3) u(6) = f(11) - f(6) u(7) = 0.5D0*f(11) - 0.5D0*f(6) - f(7) u(10) = 0.5D0*f(6) + f(7) + 1.5D0*f(11) u(11) = 0.D0 do 5 k=1,nxy3 u(k) = u(k)/ro(11) 5 continue call dumpu2(u,nx,ny) end subroutine presbd(u,ro,nx,ny,nxy3) ************************************************************* * prescribe the Neumann boundary condition. ************************************************************* integer nx,ny,nx3,nxy3,k1,k2,nb,nt,nl,nr,i,j double precision u(nxy3),ro(nxy3) double precision c1,c2,rs nx3 = nx+3 k1 = 2 k2 = nx3*(ny+2)+2 u(k1) = 1.5D0*u(k1+nx3) - 0.5D0*u(k1+nx3+1) u(k2) = 1.5D0*u(k2-nx3) - 0.5D0*u(k2-nx3+1) do 10 i=1,nx-1 k1 = k1+1 k2 = k2+1 nb = k1+nx3 nt = k2-nx3 rs = ro(k1+nx3) + ro(k1+nx3+1) c1 = ro(k1+nx3+1)/rs c2 = ro(k1+nx3)/rs u(k1) = 2.D0*u(nb) - c1*u(nb+1) - c2*u(nb-1) rs = ro(k2) + ro(k2+1) c1 = ro(k2+1)/rs c2 = ro(k2)/rs u(k2) = 2.D0*u(nt) - c1*u(nt+1) - c2*u(nt-1) 10 continue k1 = k1+1 k2 = k2+1 u(k1) = 1.5D0*u(k1+nx3) - 0.5D0*u(k1+nx3-1) u(k2) = 1.5D0*u(k2-nx3) - 0.5D0*u(k2-nx3-1) k1 = nx3+1 k2 = 2*nx3 u(k1) = 1.5D0*u(k1+1) - 0.5D0*u(k1+nx3+1) u(k2) = 1.5D0*u(k2-1) - 0.5D0*u(k2+nx3-1) do 20 j=1,ny-1 k1 = k1+nx3 k2 = k2+nx3 nl = k1+1 nr = k2-1 rs = ro(k1+nx3+1) + ro(k1+1) c1 = ro(k1+nx3+1)/rs c2 = ro(k1+1)/rs u(k1) = 2.D0*u(nl) - c1*u(nl+nx3) - c2*u(nl-nx3) rs = ro(k2) + ro(k2+nx3) c1 = ro(k2+nx3)/rs c2 = ro(k2)/rs u(k2) = 2.D0*u(nr) - c1*u(nr+nx3) - c2*u(nr-nx3) 20 continue k1 = k1+nx3 k2 = k2+nx3 u(k1) = 1.5D0*u(k1+1) - 0.5D0*u(k1+1-nx3) u(k2) = 1.5D0*u(k2-1) - 0.5D0*u(k2-1-nx3) end subroutine dumpf2(f,nx,ny) *********************************************************** c force f to be mean zero, by subtracting its average. *********************************************************** integer k,kk,nx,ny,i,j double precision f(*) double precision fsum fsum = 0.D0 k = nx+5 kk = 1 do 10 j=0,ny do 11 i=0,nx fsum = fsum + (f(k)-fsum)/dble(kk) kk = kk+1 k = k+1 11 continue k = k+2 10 continue do 20 k=1,(nx+3)*(ny+3) f(k) = f(k) - fsum 20 continue end subroutine dumpu2(u,nx,ny) ************************************************************ c force u to be mean zero, by subtracting its average. ************************************************************ integer k,kk,nx,ny,i,j double precision u(*) double precision usum,uave k(i,j) = (nx+3)*(j+1)+i+2 usum = 0.D0 do 10 j=1,ny-1 do 11 i=1,nx-1 usum = usum + u(k(i,j)) 11 continue 10 continue do 15 j=1,ny-1 usum = usum + 0.5D0*(u(k(0,j))+u(k(nx,j))) 15 continue do 16 i=1,nx-1 usum = usum + 0.5D0*(u(k(i,0))+u(k(i,ny))) 16 continue usum = usum + 0.25D0*( + u(k(0,0))+u(k(0,ny))+u(k(nx,0))+u(k(nx,ny))) uave = usum/dble(nx)/dble(ny) do 20 kk=1,(nx+3)*(ny+3) u(kk) = u(kk) - uave 20 continue end