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