program demosolve c ----------------------------------------------------------- c Basic Solver of a Poisson Equation, using MSIP c Solution Domain of Rectangular Shape c Reprogrammed for DPMS Computational Mechanics, Dec. 2006 c ----------------------------------------------------------- c implicit double precision (a-h,o-z) parameter (lx=50,ly=50,l2=lx*ly) dimension qre(l2),qrw(l2),qrn(l2),qrs(l2),qrr(l2),qrp(l2),rhs(l2) dimension ii(0:lx+1,0:ly+1),iie(l2),iiw(l2),iin(l2),iis(l2) dimension sol(0:l2),tp(0:l2) real km dimension am(0:l2),bm(0:l2),cm(0:l2),dm(0:l2),em(0:l2), 1 fm(0:l2),gm(0:l2),hm(0:l2),km(0:l2) c c Data Input c ----------- write(*,*)' ENTER GRID DIMENSIONS nx ny ' read(*,*)nx,ny write(*,*)' DIMENSIONS : [ ',nx,' x ',ny,' ] ' write(*,*) ' ENTER Number of Iterations & PSI (or alfa) value' read(*,*)itgmr,psi c ismax=nx*ny if(ismax.gt.l2) stop ' INSUFFICIENT DIMENSIONING ' c c Topology and Connectivities c --------------------------- call string(lx,ly,l2,nx,ny,ismax,ii,iie,iiw,iin,iis) c c Linear Problem Definition c ------------------------- c-The problem is linear, so MATRI is called once! call matri(lx,ly,l2,nx,ny,ismax,ii,iie,iiw,iin,iis, 1 qre,qrw,qrn,qrs,qrr,qrp,rhs) c c Factorization (MSIP5 or MSIP9) c ------------------------------ call msip9(lx,ly,l2,nx,ny,ismax,ii,iie,iiw,iin,iis, 1 qre,qrw,qrn,qrs,qrr,qrp,rhs,psi, 2 am,bm,cm,dm,em,fm,gm,hm,km) c c Iterative Solver c ---------------- call solver(lx,ly,l2,nx,ny,ismax,ii,iie,iiw,iin,iis, 1 qre,qrw,qrn,qrs,qrr,qrp,rhs, 2 am,bm,cm,dm,em,fm,gm,hm,km,sol,itgmr,tp) c c Results Printout c ---------------- open(1,file='res') write(1,*)' Solution obtained using PSI= ',psi write(1,*)' ------------------ ' write(1,*)nx,ny do 80 ix=1,nx do 80 iy=1,ny l=ii(ix,iy) xl=dfloat(ix) yl=dfloat(iy) solu=dfloat(ix)**2+dfloat(iy)**2 ! exact solution write(1,6)xl,yl,sol(l),solu 6 format(4(1x,f16.7)) 80 continue close(1) c stop end c c c*****************************************************************c subroutine string(lx,ly,l2,nx,ny,ismax,ii,iie,iiw,iin,iis) c*****************************************************************c implicit double precision (a-h,o-z) dimension ii(0:lx+1,0:ly+1),iie(l2),iiw(l2),iin(l2),iis(l2) c c-Determines the topology-connectivity (in case of complex shapes) c-So: ii(ix,iy)=string number of node (ix,iy) c- and for internal nodes: c- iie=east(ix+1),iiw=west(ix-1),iin=north(iy+1),iis=south(iy-1) c- whereas for boundary nodes iie or iiw or iis or iin c- are set to zero to denote that this node does not exist c do 3 ix=0,nx+1 do 3 iy=0,ny+1 ii(ix,iy)=0 3 continue c ismax=0 ! 1D counter of nodes do 10 ix=1,nx do 10 iy=1,ny ismax=ismax+1 ii(ix,iy)=ismax 10 continue write(*,*)' ismax = ',ismax c do 15 ix=1,nx do 15 iy=1,ny L=ii(ix,iy) if (ix.lt.nx) iie(L)=ii(ix+1,iy) if (ix.gt.1) iiw(L)=ii(ix-1,iy) if (iy.lt.ny) iin(L)=ii(ix,iy+1) if (iy.gt.1) iis(L)=ii(ix,iy-1) 15 continue c return end c c c ************************************************************ subroutine matri(lx,ly,l2,nx,ny,ismax,ii,iie,iiw,iin,iis, 1 qre,qrw,qrn,qrs,qrr,qrp,rhs) c ************************************************************ implicit double precision (a-h,o-z) dimension qre(l2),qrw(l2),qrn(l2),qrs(l2),qrr(l2),qrp(l2),rhs(l2) dimension ii(0:lx+1,0:ly+1),iie(l2),iiw(l2),iin(l2),iis(l2) c c-Fills in the coefficient matrix and the rhs array for c-two demo problems with closed form solutions. c-Be careful, the equation to be solved reads as follows: c- QRE(L)*PHI(IIE) + QRW(L)*PHI(IIW) + c- + QRN(L)*PHI(IIN) + QRS(L)*PHI(IIS) - c- - QRP(L)*PHI( L ) + RHS(L) = 0 c-(Please note the signs of the last two terms) c- c-The demo problem is: c- LAPLACE(x^2+y^2)=4 solved on a simple orthogonal grid with c- Delta_x=Delta_y=1, with Dirichlet boundary conditions c do 1 ix=1,nx do 1 iy=1,ny L=ii(ix,iy) if(ix.eq.1.or.ix.eq.nx.or.iy.eq.1.or.iy.eq.ny) then qrp(L) = 1.d0 qre(L) = 0.d0 qrw(L) = 0.d0 qrn(L) = 0.d0 qrs(L) = 0.d0 rhs(L) = dfloat(ix)**2+dfloat(iy)**2 else qre(L) = 1.d0 qrw(L) = 1.d0 qrn(L) = 1.d0 qrs(L) = 1.d0 qrp(L)=qre(L)+qrw(L)+qrn(L)+qrs(L) rhs(L) = -4.d0 endif 1 continue c return end c c c ************************************************************ subroutine msip5(lx,ly,l2,nx,ny,ismax,ii,iie,iiw,iin,iis, 1 qre,qrw,qrn,qrs,qrr,qrp,rhs,psi, 2 am,bm,cm,dm,em,fm,gm,hm,km) c ************************************************************ implicit double precision (a-h,o-z) dimension qre(l2),qrw(l2),qrn(l2),qrs(l2),qrr(l2),qrp(l2),rhs(l2) dimension ii(0:lx+1,0:ly+1),iie(l2),iiw(l2),iin(l2),iis(l2) real km,kc dimension am(0:l2),bm(0:l2),cm(0:l2),dm(0:l2),em(0:l2), 1 fm(0:l2),gm(0:l2),hm(0:l2),km(0:l2) c c- (iy) C -- F -- K c- | B --[E]-- H ---(ix) c- | A -- D -- G c c c Initialization c -------------- do 4 ix=0,nx+1 do 4 iy=0,ny+1 L=ii(ix,iy) am(L)=0.d0 bm(L)=0.d0 cm(L)=0.d0 dm(L)=0.d0 em(L)=0.d0 fm(L)=0.d0 gm(L)=0.d0 hm(L)=0.d0 km(L)=0.d0 4 continue c c MSIP, stencil 5 nodes c --------------------- do 1 ix=1,nx do 1 iy=1,ny L=ii(ix,iy) ie=iie(L) iw=iiw(L) in=iin(L) is=iis(L) c-Capital Letters (according to the distributes notes) ec=-qrp(L) hc= qre(L) bc= qrw(L) fc= qrn(L) dc= qrs(L) kc= 0.d0 cc= 0.d0 gc= 0.d0 ac= 0.d0 c-Small Letters (according to the distributes notes) bm(L) = bc / (1.d0+psi*fm(iw)) dm(L) = dc / (1.d0+psi*hm(is)) em(L) = ec - bm(L)*hm(iw) - dm(L)*fm(is) 1 + psi * ( bm(L)*fm(iw) + dm(L)*hm(is) ) if (dabs(em(L)).lt..1e-10) stop ' Stop in Factor' em(L) = 1.d0/em(L) ! attention keeps the INVERSE of EPSILON! fm(L) = ( fc - psi*bm(L)*fm(iw) ) * em(L) hm(L) = ( hc - psi*dm(L)*hm(is) ) * em(L) 1 continue c return end c c c ************************************************************ subroutine msip9(lx,ly,l2,nx,ny,ismax,ii,iie,iiw,iin,iis, 1 qre,qrw,qrn,qrs,qrr,qrp,rhs,psi, 2 am,bm,cm,dm,em,fm,gm,hm,km) c ************************************************************ implicit double precision (a-h,o-z) dimension qre(l2),qrw(l2),qrn(l2),qrs(l2),qrr(l2),qrp(l2),rhs(l2) dimension ii(0:lx+1,0:ly+1),iie(l2),iiw(l2),iin(l2),iis(l2) real km,kc dimension am(0:l2),bm(0:l2),cm(0:l2),dm(0:l2),em(0:l2), 1 fm(0:l2),gm(0:l2),hm(0:l2),km(0:l2) c c- (iy) C -- F -- K c- | B --[E]-- H ---(ix) c- | A -- D -- G c c c Initialization c -------------- do 4 ix=0,nx+1 do 4 iy=0,ny+1 L=ii(ix,iy) am(L)=0.d0 bm(L)=0.d0 cm(L)=0.d0 dm(L)=0.d0 em(L)=0.d0 fm(L)=0.d0 gm(L)=0.d0 hm(L)=0.d0 km(L)=0.d0 4 continue c c MSIP, stencil 5 nodes c --------------------- do 1 ix=1,nx do 1 iy=1,ny L=ii(ix,iy) ie=iie(L) iw=iiw(L) in=iin(L) is=iis(L) if(in.ne.0.and.ie.ne.0) then ine=iin(ie) else ine=0 endif if(is.ne.0.and.ie.ne.0) then ise=iis(ie) else ise=0 endif if(in.ne.0.and.iw.ne.0) then inw=iin(iw) else inw=0 endif if(is.ne.0.and.iw.ne.0) then isw=iis(iw) else isw=0 endif c-Capital Letters (according to the distributes notes) ec=-qrp(L) hc= qre(L) bc= qrw(L) fc= qrn(L) dc= qrs(L) kc= 0.d0 cc= 0.d0 gc= 0.d0 ac= 0.d0 c-Small Letters (according to the distributes notes) am(L)=ac bm(L)=(bc-psi*cc*fm(inw)-am(L)*fm(isw))/(1.d0-psi*fm(iw)*fm(inw)) cm(L)=cc-bm(L)*fm(iw) dm(L)=(dc-am(L)*(2.*psi*gm(isw)+hm(isw))-bm(L)*gm(iw)) 1 / (1.d0+2.*psi*gm(is)) em(L)=ec+am(L)*(psi*gm(isw)-km(isw))-bm(L)*hm(iw)+ 1 cm(L)*(2.*psi*fm(inw)-gm(inw)+psi*km(inw))+ 2 dm(L)*(2.*psi*gm(is)-fm(is)) if (dabs(em(L)).lt..1e-10) stop ' Stop in Factor' em(L) = 1.d0/em(L) ! attention keeps the INVERSE of EPSILON! fm(L)=( fc-bm(L)*km(iw)-cm(L)* 1 (hm(inw)+2.*psi*fm(inw)+2.*psi*km(inw))) * em(L) gm(L)=(gc-dm(L)*hm(is)) * em(L) hm(L)=(hc-dm(L)*(psi*gm(is)+km(is))) * em(L) km(L)=kc * em(L) 1 continue c return end c c c ************************************************************ subroutine solver(lx,ly,l2,nx,ny,ismax,ii,iie,iiw,iin,iis, 1 qre,qrw,qrn,qrs,qrr,qrp,rhs, 2 am,bm,cm,dm,em,fm,gm,hm,km,sol,itgmr,tp) c ************************************************************ implicit double precision (a-h,o-z) dimension qre(l2),qrw(l2),qrn(l2),qrs(l2),qrr(l2),qrp(l2),rhs(l2) dimension ii(0:lx+1,0:ly+1),iie(l2),iiw(l2),iin(l2),iis(l2) real km dimension am(0:l2),bm(0:l2),cm(0:l2),dm(0:l2),em(0:l2), 1 fm(0:l2),gm(0:l2),hm(0:l2),km(0:l2) dimension sol(0:l2),tp(0:l2) c itmmn=1 open(2,file='conv.dat') c 3 do 1 iter=itmmn,itgmr c c Delta Formulation c ----------------- resit=0.d0 do L=1,ismax ie=iie(L) iw=iiw(L) in=iin(L) is=iis(L) qrr(L) = - ( qre(L)*sol(ie)+qrw(L)*sol(iw)+ 1 qrn(L)*sol(in)+qrs(L)*sol(is) 1 -qrp(L)*sol(L)+rhs(L) ) resit=resit+qrr(L)*qrr(L) enddo resit=dsqrt(resit)/dfloat(ismax) c write(*,7)iter-1,resit write(2,8)iter-1,resit 7 format(2x,' ***** ITERATION : ',i5,' , RESIDUAL : ',e14.7) 8 format(2x,i5,e14.7) c call backfron(lx,ly,l2,nx,ny,ismax,ii,iie,iiw,iin,iis, 1 qrr,am,bm,cm,dm,em,fm,gm,hm,km,sol,tp) c do L=1,ismax sol(L)=sol(L)+tp(L) enddo c 1 continue c write(*,*)'ENTER 0 TO STOP OR ITMAX TO CONTINUE' read(*,*)itmix if(itmix.ne.0)then itmmn=itgmr+1 itgmr=itgmr+itmix go to 3 endif c close(2) return end c c c ************************************************************ subroutine backfron(lx,ly,l2,nx,ny,ismax,ii,iie,iiw,iin,iis, 1 qrr,am,bm,cm,dm,em,fm,gm,hm,km,sol,tp) c ************************************************************ implicit double precision (a-h,o-z) dimension qrr(l2),tp(0:l2),aux(0:l2),sol(0:l2) dimension ii(0:lx+1,0:ly+1),iie(l2),iiw(l2),iin(l2),iis(l2) real km dimension am(0:l2),bm(0:l2),cm(0:l2),dm(0:l2),em(0:l2), 1 fm(0:l2),gm(0:l2),hm(0:l2),km(0:l2) c do L=0,ismax aux(L) = 0.d0 tp(L) = 0.d0 enddo c c Front Substitution c ------------------ do 10 ix=1,nx do 10 iy=1,ny L=ii(ix,iy) iw=iiw(L) is=iis(L) in=iin(L) if(in.ne.0.and.iw.ne.0) then inw=iin(iw) else inw=0 endif if(is.ne.0.and.iw.ne.0) then isw=iis(iw) else isw=0 endif comw = 0.d0 comsw = 0.d0 comnw = 0.d0 coms = 0.d0 if(isw.ne.0) comsw = am(L)*aux(isw) if(iw.ne.0) comw = bm(L)*aux(iw) if(inw.ne.0) comnw = cm(L)*aux(inw) if(is.ne.0) coms = dm(L)*aux(is) aux(L)=(qrr(L)-comsw-comw-comnw-coms) * em(L) 10 continue c c Back Substitution c ----------------- do 20 ix=nx,1,-1 do 20 iy=ny,1,-1 L=ii(ix,iy) in=iin(L) ie=iie(L) if(in.ne.0.and.ie.ne.0) then ine=iin(ie) else ine=0 endif is=iis(L) if(is.ne.0.and.ie.ne.0) then ise=iis(ie) else ise=0 endif comn = 0.d0 comse = 0.d0 come = 0.d0 comne = 0.d0 if(in.ne.0) comn = fm(L)*tp(in) if(ise.ne.0) comse = gm(L)*tp(ise) if(ie.ne.0) come = hm(L)*tp(ie) if(ine.ne.0) comne = km(L)*tp(ine) tp(L)=aux(L)-comn-comse-come-comne 20 continue c return end c c