subroutine poisson(f,u,res,bl,br,bb,bt,nx,ny,mxy) ***************************************************************** * * MULTIGRID PACKAGE TO SOLVE THE 5-POINT CONSTANT COEFFECIENT * POISSON EQUATION: * * div grad u = f * * 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),res(mxy) double precision bl(ny),br(ny),bb(nx),bt(nx) double precision el2,h,hinv,hh,eps common/meshsz/h,hinv parameter (kit=2,nv=8) parameter (eps=1.D-10) call modify(f,bl,br,bb,bt,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) 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),kit) call resid1(nnx,nny,nxy2,hh,f(kf),u(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),kit) 26 continue call resid1(nx,ny,nxy2,h,f,u,res,el2) c print *,'cycle ',kv,' el2= ',el2 if (el2 .le. eps) goto 6 call dumpnu(u,nx,ny) 5 continue 6 print *,'poisson solver el2',el2 101 format(2(f12.6,2x)) end subroutine modify(f,bl,br,bb,bt,nx,ny) **************************************************************** * incorporate the Neumann boundary data into the right-hand- * side f **************************************************************** integer nx,ny,i,j double precision f(0:nx+1,0:ny+1) double precision bl(ny),br(ny),bb(nx),bt(ny) double precision h,hinv common/meshsz/h,hinv do 10 j=1,ny do 10 i=1,nx f(i,j) = f(i,j)*h*h 10 continue do 20 j=1,ny f(1,j) = f(1,j) - bl(j)*h f(nx,j) = f(nx,j) - br(j)*h 20 continue do 30 i=1,nx f(i,1) = f(i,1) - bb(i)*h f(i,ny) = f(i,ny) - bt(i)*h 30 continue call dumpnu(f,nx,ny) end subroutine gsneu1(nx,ny,nxy2,u,f,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) 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 u(k) = 0.25D0*(u(k+1)+u(k-1)+u(k+nx2)+u(k-nx2)-f(k)) 30 continue do 40 j=2,ny,2 do 40 i=2,nx,2 k = nx2*j+i+1 u(k) = 0.25D0*(u(k+1)+u(k-1)+u(k+nx2)+u(k-nx2)-f(k)) 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 u(k) = 0.25D0*(u(k+1)+u(k-1)+u(k+nx2)+u(k-nx2)-f(k)) 50 continue do 60 j=2,ny,2 do 60 i=1,nx-1,2 k = nx2*j+i+1 u(k) = 0.25D0*(u(k+1)+u(k-1)+u(k+nx2)+u(k-nx2)-f(k)) 60 continue 5 continue end subroutine resid1(nx,ny,nxy2,h,f,u,res,el2) ****************************************************************** * compute the residue ******************************************************* integer nx,ny,i,j,k,nx2,nxy2 double precision f(nxy2),res(nxy2),u(nxy2) double precision el2,h nx2 = nx+2 call presb1(u,nx,ny,nxy2) do 10 j=1,ny do 12 i=1,nx k = nx2*j+i+1 res(k) = f(k) + 4.D0*u(k) - + u(k+1)-u(k-1)-u(k+nx2)-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 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 setzro(u,nxy) integer k,nxy double precision u(nxy) do 10 k=1,nxy u(k) = 0.D0 10 continue end