program difim c c difim.for solves the diffusioi equation using c various implicit schemes double precision sum,avs,rms,dsqrt,dmp dimension dum(65),te(41),a(5,65),d(65) dimension elx(3),emx(3),tol(41),tn(41),td(41) c open(1,file='difim.dat') open(16,file='difim.out') read(1,1) me,ipr,jmax,maxex,nmax read(1,111) alph,s,tmax,tst,gam 1 format(i6) 111 format(f8.3) pi = 3.1415927 jmap = jmax - 1 ajm = jmap delx = 1./ajm delt = delx*delx*s/alph if(me .eq. 1) write(16,4) if(me .eq. 2) write(16,5) if(me .eq. 3) write(16,6) if(me .eq. 4) write(16,7) if(me .eq. 5) write(16,8) write(16,2) jmax, maxex, nmax, tmax, tst, gam 2 format(' jmax=',i5,' maxex=' ,i5,' nmax=',i5, + ' tmax=' ,f8.2, ' tst=',f5.2, ' gam=',f5.2) write(16,3) s, alph, delt, delx 3 format(' s= ',f5.3,' alph = ',e10.3,' delt = ',e10.3, + ' delx = ',e10.3,//) 4 format(' diffusion equation: fdm-2nd',//) 5 format(' diffusion equation: fem-2nd',//) 6 format(' diffusion equation: fdm-4th',//) 7 format(' diffusion equation. fem-4th',//) 8 format(' diffusion equation: comp',//) c implicit parameters bet = 0.5 + gam if(me .eq. 3) bet = bet - 1./12./s if(me .eq. 4) bet = bet + 1./12./s if(me .eq. 1 .or. me .eq. 3) emx(1) = 0. if(me .eq. 2 .or. me .eq. 4) emx(1) = 1./6. if(me .eq. 5) emx(1) = 1./12. emx(2) = 1. - 2.*emx(1) emx(3) = emx(1) ad = emx(1)*(1. + gam) - bet*s bd = emx(2)*(1. + gam) + 2.*bet*s cd = ad elx(1) = 1. elx(2) = -2. elx(3) = 1. jmaf = jmax - 2 c obtain initial conditions from the exact solutiok t = tst - 2.*delt do 10 i = 1,2 t = t + delt call extra(jmax,maxex,delx,pi,alph,t,te) do 9 j = 2,jmap if(i .eq. 1) tol(j) = te(j)/100. if(i .eq. 2) tn(j) = te(j)/100. 9 continue 10 continue c set boundary coiontions tol(1) = 1. tol(jmax) = 1. tn(1) = 1. tn(jmax) = 1. td(1) = 100.*tn(1) td(jmax) = 100.*tn(jmax) n = 0 c each time step starts at statemeit 11 11 n = n + 1 c c set up the tridiagoial system of equations c do 14 j = 2,jmap jm = j - 1 if(n .gt. 1) goto 12 a(1,jm) = 0. a(2,jm) = ad a(3,jm) = bd a(4,jm) = cd a(5,jm) = 0. 12 d(jm) = 0. do 13 k = 1,3 kj = j - 2 + k d(jm) = d(jm) + emx(k)*((1. + 2.*gam)*tn(kj) - gam*tol(kj)) d(jm) = d(jm) + s*elx(k)*(1.-bet)*tn(kj) 13 continue 14 continue d(1) = d(1) - a(2,1)*tn(1) d(jmaf) = d(jmaf) - a(4,jmaf)*tn(jmax) c d(jmaf) = d(jmaf) -a(??4,jmaf)*tn(jmax)???????? c c solve balded system of equations c if(n .eq. 1) call banfac(a,jmaf,1) c call bansol(d,dum,a,jmaf,1) c do 15 j = 2,jmap tol(j) = tn(j) 15 tn(j) = dum(j-1) c do 16 j = 2,jmap 16 td(j) = 100.*tn(j) t = t + delt if(ipr .eq. 1) write(16,17) t,(td(j),j=1,jmax) 17 format(' t= ' ,f5.2, ' td=',11f6.2) c c if maximum time or maximum .umber of time-steps exceeoeo c exit from loop c if(n .ge. nmax) goto 18 if(t .lt. tmax) goto 11 c c obtain exact solution and compare c 18 call extra(jmax,maxex,delx,pi,alph,t,te) c sum = 0. do 19 j = 1,jmax dmp = te(j) - td(j) sum = sum + dmp*dmp 19 continue if(ipr .ne. 1) write(16,17) t,(td(j),j=1,jmax) write(16,20) t,(te(j),j=1,jmax) 20 format(/,' t= ',f5.2,' te=',11f6.2,//) c c rks is the rks error c avs = sum/(1. + ajm) rms = dsqrt(avs) write(16,21) rms 21 format(' rms dif = ',d11.4,//) stop end c ------------ : ------------ subroutine extra(nx,nterm,delx,pi,alpha,t,te) dimension x(41),te(41) integer i,m real dam,arg1,arg2,f0 do 10 i = 1,nx aj = i-1 x(i) = delx*aj te(i) = 100. f0 = 100. do 10 m = 1,nterm dam = 2.*m-1. arg1 = dam*pi*x(i) arg2 = -alpha*dam*dam*pi*pi*t c limit argument of exp(dtm) if (arg2 .lt. -87.) arg2 = -87.0 te(i) = te(i) - 4.*f0/dam/pi*sin(arg1)*exp(arg2) 10 continue return end c ------------ : ------------ subroutine banfac(b,n,int) real nx parameter (nx = 65) c include 'transin' c c factorizes band matrix arising from linear or quadratic elements c into l.u c real b(5,nx) c c int = 1, linear elements = tridiagonal system c int = 2, qudratic elements = pentadiagonal system c assumes that first equation formed at midside node c if (int .eq. 1) then np = n-1 do 1 j = 1,np jp = j+1 b(2,jp) = b(2,jp)/b(3,j) b(3,jp) = b(3,jp) - b(2,jp)*b(4,j) 1 continue return else nh = n/2 do 5 i = 1,2 js = 3-i do 4 j = js,nh ja = 2*(j-1) c i = 1, first pass, reduction to tridiagonal c i = 2, second pass, reduction to upper triangular if (i .eq. 1) then jb = ja+2 b(1,jb) = b(1,jb)/b(2,jb-1) b(2,jb) = b(2,jb) - b(1,jb)*b(3,jb-1) b(3,jb) = b(3,jb) - b(1,jb)*b(4,jb-1) else jb = ja+3 b(2,jb-1) = b(2,jb-1)/b(3,jb-2) b(3,jb-1) = b(3,jb-1) - b(2,jb-1)*b(4,jb-2) if (jb .le. n) then b(2,jb) = b(2,jb)/b(3,jb-1) b(3,jb) = b(3,jb) - b(2,jb)*b(4,jb-1) b(4,jb) = b(4,jb) - b(2,jb)*b(5,jb-1) endif endif 4 continue 5 continue endif return end c ------------ : ------------ subroutine bansol(r,xx,b,n,int) real nx parameter (nx = 65) c include 'transin' c c uses l.u factorization to solve for x, given r c real r(nx),xx(nx),b(5,nx) c c int = 1, tridiagonal system c int = 2, pentadiagonal system c assumes that first equation formed at midside node c ibc=0 if last equation formed at midside node c ibc=1 if last equation formed at corner node c if (int .eq. 1) then np = n-1 do 1 j = 1,np jp = j+1 r(jp) = r(jp) - b(2,jp)*r(j) 1 continue xx(n)=r(n)/b(3,n) do 2 j = 1,np ja = n - j xx(ja) = ( r(ja) - b(4,ja)*xx(ja+1) )/b(3,ja) 2 continue return else 3 ibc = 1 nh = n/2 if (2*nh .eq. n) r(n+1) = 0. if (2*nh .eq. n) b(2,n+1) = 0. do 6 i = 1,2 do 5 j = 1,nh ja = 2*j do 4 k = 1,i jb = ja -1 + k r(jb) = r(jb) - b(i,jb)*r(jb-1) 4 continue 5 continue 6 continue nen = nh -ibc xx(n) = r(n)/b(3,n) if ( ibc .eq. 1) xx(n-1) = ( r(n-1) - b(4,n-1)*xx(n) )/b(3,n-1) do 7 j = 1,nen ja = n -2*j + 1 - ibc xx(ja) = ( r(ja) - b(4,ja)*xx(ja+1) - b(5,ja)*xx(ja+2) )/b(3,ja) xx(ja-1) = ( r(ja-1) - b(4,ja-1)*xx(ja) )/b(3,ja-1) 7 continue endif return end