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