program demosolve c ----------------------------------------------------------- c Basic Solver of a simple Poisson Equation, using MSIP c For the JPC Computational Mechanics, Aug. 2016 c ----------------------------------------------------------- c implicit double precision (a-h,o-z) parameter (lx=50,ly=50,l2=lx*ly) dimension ii(0:lx+1,0:ly+1) dimension sol(0:l2),dsol(0:l2),x(0:l2),y(0:l2),rhs(l2),qrr(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 if(nx.gt.lx.or.ny.gt.ly) stop ' Increase LX or/and LY ' write(*,*)' DIMENSIONS : [ ',nx,' x ',ny,' ] ' write(*,*) ' ENTER Number of Iterations & PSI (or alfa) value' read(*,*)itgmr,psi c c String Numbering c ---------------- call string(lx,ly,l2,nx,ny,ismax,ii) c c Grid Coordinates c ---------------- do ix=1,nx do iy=1,ny L=ii(ix,iy) x(L)=dfloat(ix) y(L)=dfloat(iy) sol(L)=0.d0 enddo enddo c c Computation of the RHS term & Dirichlet Conditions c -------------------------------------------------- call rhsDiri(lx,ly,l2,nx,ny,ismax,ii,x,y,rhs,sol) c c Factorization (MSIP5 or MSIP9) c ------------------------------ call msip9(lx,ly,l2,nx,ny,ismax,ii,x,y,psi,sol, 2 am,bm,cm,dm,em,fm,gm,hm,km) c c Iterative Solver (Back-Front Substitution) c ------------------------------------------ call solver(lx,ly,l2,nx,ny,ismax,ii,x,y,qrr,rhs, 2 am,bm,cm,dm,em,fm,gm,hm,km,sol,itgmr,dsol) c c Results Printout c ---------------- open(1,file='res') write(1,*)' Solution obtained using PSI= ',psi write(1,*)' ------------------ ' write(1,*)nx,ny do ix=1,nx do iy=1,ny L=ii(ix,iy) xl=x(L) yl=y(L) solu=x(L)**2+y(L)**2 ! exact solution write(1,'(4(1x,f22.12))')xl,yl,sol(l),solu enddo enddo close(1) c stop end c c c*****************************************************************c subroutine string(lx,ly,l2,nx,ny,ismax,ii) c*****************************************************************c implicit double precision (a-h,o-z) dimension ii(0:lx+1,0:ly+1) c do ix=0,nx+1 do iy=0,ny+1 ii(ix,iy)=0 enddo enddo c ismax=0 ! 1D counter of nodes do ix=1,nx do iy=1,ny ismax=ismax+1 ii(ix,iy)=ismax enddo enddo write(*,*)' ismax = ',ismax c return end c c c ************************************************************ subroutine rhsDiri(lx,ly,l2,nx,ny,ismax,ii,x,y,rhs,sol) c ************************************************************ implicit double precision (a-h,o-z) dimension rhs(l2),x(0:l2),y(0:l2),sol(0:l2) dimension ii(0:lx+1,0:ly+1) c do ix=1,nx do iy=1,ny L=ii(ix,iy) if(ix.eq.1.or.ix.eq.nx.or.iy.eq.1.or.iy.eq.ny) then rhs(L) = x(L)**2+y(L)**2 sol(L) = rhs(L) else rhs(L) = -4.d0 endif enddo enddo c return end c c c ************************************************************ subroutine msip5(lx,ly,l2,nx,ny,ismax,ii,x,y,psi,sol, 2 am,bm,cm,dm,em,fm,gm,hm,km) c ************************************************************ implicit double precision (a-h,o-z) dimension ii(0:lx+1,0:ly+1) real km,kc dimension x(0:l2),y(0:l2),sol(0:l2) 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 L=0,ismax bm(L)=0.d0 dm(L)=0.d0 em(L)=0.d0 fm(L)=0.d0 hm(L)=0.d0 enddo c c MSIP, stencil 9 nodes c --------------------- do ix=1,nx do iy=1,ny L=ii(ix,iy) ie=ii(ix+1,iy) iw=ii(ix-1,iy) in=ii(ix,iy+1) is=ii(ix,iy-1) c-Capital Letters (see course notes) if(ix.eq.1.or.ix.eq.nx.or.iy.eq.1.or.iy.eq.ny) then ec= 1.d0 hc= 0.d0 bc= 0.d0 fc= 0.d0 dc= 0.d0 else ec=-4.d0 hc= 1.d0 bc= 1.d0 fc= 1.d0 dc= 1.d0 endif c-Small Letters (see course notes), elements of L and U 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) enddo enddo c return end c c c ************************************************************ subroutine msip9(lx,ly,l2,nx,ny,ismax,ii,x,y,psi,sol, 2 am,bm,cm,dm,em,fm,gm,hm,km) c ************************************************************ implicit double precision (a-h,o-z) dimension ii(0:lx+1,0:ly+1) real km,kc dimension x(0:l2),y(0:l2),sol(0:l2) 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 L=0,ismax 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 enddo c c MSIP, stencil 9 nodes c --------------------- do ix=1,nx do iy=1,ny L=ii(ix,iy) ie=ii(ix+1,iy) iw=ii(ix-1,iy) in=ii(ix,iy+1) is=ii(ix,iy-1) ine=ii(ix+1,iy+1) ise=ii(ix+1,iy-1) inw=ii(ix-1,iy+1) isw=ii(ix-1,iy-1) c-Capital Letters (see course notes) if(ix.eq.1.or.ix.eq.nx.or.iy.eq.1.or.iy.eq.ny) then ec= 1.d0 hc= 0.d0 bc= 0.d0 fc= 0.d0 dc= 0.d0 kc= 0.d0 cc= 0.d0 gc= 0.d0 ac= 0.d0 else ec=-4.d0 hc= 1.d0 bc= 1.d0 fc= 1.d0 dc= 1.d0 kc= 0.d0 cc= 0.d0 gc= 0.d0 ac= 0.d0 endif c-Small Letters (see course notes), elements of L and U 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) enddo enddo c return end c c c ************************************************************ subroutine solver(lx,ly,l2,nx,ny,ismax,ii,x,y,qrr,rhs, 2 am,bm,cm,dm,em,fm,gm,hm,km,sol,itgmr,dsol) c ************************************************************ implicit double precision (a-h,o-z) dimension qrr(l2),rhs(l2) dimension ii(0:lx+1,0:ly+1) 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),dsol(0:l2),x(0:l2),y(0:l2) c itmmn=1 open(2,file='conv.dat') c 3 do iter=itmmn,itgmr c c Delta Formulation - Compute Residual c ------------------------------------ resit=0.d0 do ix=1,nx do iy=1,ny L=ii(ix,iy) ie=ii(ix+1,iy) iw=ii(ix-1,iy) in=ii(ix,iy+1) is=ii(ix,iy-1) ine=ii(ix+1,iy+1) ise=ii(ix+1,iy-1) inw=ii(ix-1,iy+1) isw=ii(ix-1,iy-1) if(ix.eq.1.or.ix.eq.nx.or.iy.eq.1.or.iy.eq.ny) then qrr(L) = sol(L)-rhs(L) else qrr(L) = -(sol(ie)+sol(iw)+sol(in)+sol(is)-4.d0*sol(L)+rhs(L)) endif resit=resit+qrr(L)*qrr(L) enddo 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 c Back-Front Substitution c ----------------------- call backfron(lx,ly,l2,nx,ny,ismax,ii, 1 qrr,am,bm,cm,dm,em,fm,gm,hm,km,sol,dsol) c do L=1,ismax sol(L)=sol(L)+dsol(L) enddo c enddo 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, 1 qrr,am,bm,cm,dm,em,fm,gm,hm,km,sol,dsol) c ************************************************************ implicit double precision (a-h,o-z) dimension qrr(l2),dsol(0:l2),aux(0:l2),sol(0:l2) dimension ii(0:lx+1,0:ly+1) 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 dsol(L) = 0.d0 enddo c c Front Substitution c ------------------ do ix=1,nx do iy=1,ny L=ii(ix,iy) iw=ii(ix-1,iy) is=ii(ix,iy-1) inw=ii(ix-1,iy+1) isw=ii(ix-1,iy-1) comsw = am(L)*aux(isw) comw = bm(L)*aux(iw) comnw = cm(L)*aux(inw) coms = dm(L)*aux(is) aux(L)=(qrr(L)-comsw-comw-comnw-coms) * em(L) enddo enddo c c Back Substitution c ----------------- do ix=nx,1,-1 do iy=ny,1,-1 L=ii(ix,iy) ie=ii(ix+1,iy) in=ii(ix,iy+1) ine=ii(ix+1,iy+1) ise=ii(ix+1,iy-1) comn = fm(L)*dsol(in) comse = gm(L)*dsol(ise) come = hm(L)*dsol(ie) comne = km(L)*dsol(ine) dsol(L)=aux(L)-comn-comse-come-comne enddo enddo c return end c c