c fishpk37 from portlib 12/30/83 c c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * c * * c * f i s h p a k * c * * c * * c * a package of fortran subprograms for the solution of * c * * c * separable elliptic partial differential equations * c * * c * (version 3.1 , october 1980) * c * * c * by * c * * c * john adams, paul swarztrauber and roland sweet * c * * c * of * c * * c * the national center for atmospheric research * c * * c * boulder, colorado (80307) u.s.a. * c * * c * which is sponsored by * c * * c * the national science foundation * c * * c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * c c c an example showing the use of sepeli to solve the separable c elliptic partial differential equation . . . c (x+1)**2*uxx+2*(x+1)*ux+exp(y)*uyy-(x+y)*u = g(x,y) on c 0.le.x.le.1, 0.le.y.le.1 with specified boundary conditions c at y=0,1 and mixed boundary conditions of the form c ux(0,y)+u(0,y), ux(1,y)+u(1,y) at x=0,1. c the approximation is generated on a uniform 33 by 33 grid. c the exact solution u(x,y)=(x*y)**3+1 is used to set the c right hand side, boundary conditions, and compute second and c fourth order discretization error c the exact work space length required is 1118 words. c this was determined by a previous call to sepeli and print c out of w(1). c dimension usol(33,33),grhs(33,33),bda(33) ,bdb(33) , 1 w(1118) c c declare coefficient subroutines external c external cofx ,cofy c c define arithmetic functions giving exact solution c ue(s,t) = (s*t)**3+1.0 uxe(s,t) = 3.0*s**2*t**3 uxxe(s,t) = 6.0*s*t**3 uye(s,t) = 3.0*s**3*t**2 uyye(s,t) = 6.0*s**3*t c c set limits on region c a = 0.0 b = 1.0 c = 0.0 d = 1.0 c c set grid size c m = 32 n = 32 dlx = (b-a)/float(m) dly = (d-c)/float(n) nx = m+1 ny = n+1 do 102 i=1,nx x = a+float(i-1)*dlx c c set specified boundary conditions at y=c,d c usol(i,1) = ue(x,c) usol(i,ny) = ue(x,d) call cofx (x,af,bf,cf) do 101 j=1,ny y = c+float(j-1)*dly call cofy (y,df,ef,ff) c c set right hand side c grhs(i,j) = af*uxxe(x,y)+bf*uxe(x,y)+cf*ue(x,y)+ 1 df*uyye(x,y)+ef*uye(x,y)+ff*ue(x,y) 101 continue 102 continue c c set mixed boundary conditions at x=a,b c alpha = 1.0 beta = 1.0 do 103 j=1,ny y = c+float(j-1)*dly bda(j) = uxe(a,y)+alpha*ue(a,y) bdb(j) = uxe(b,y)+beta*ue(b,y) 103 continue c c set boundary swithces c mbdcnd = 3 nbdcnd = 1 c c set first dimension of usol,grhs and work space length c idmn = 33 w(1) = 1118. c c set work space length in first word c set initial call parameter to zero c intl = 0 c c obtain second order approximation c iorder = 2 call sepeli (intl,iorder,a,b,m,mbdcnd,bda,alpha,bdb,beta,c,d,n, 1 nbdcnd,dum,dum,dum,dum,cofx,cofy,grhs,usol,idmn,w, 2 pertrb,ierror) err = 0.0 do 105 i=1,nx x = a+float(i-1)*dlx do 104 j=1,ny y = c+float(j-1)*dly err = amax1(err,abs((usol(i,j)-ue(x,y))/ue(x,y))) 104 continue 105 continue err2 = err c c obtain fourth order approximation c iorder = 4 c c non-initial call c intl = 1 call sepeli (intl,iorder,a,b,m,mbdcnd,bda,alpha,bdb,beta,c,d,n, 1 nbdcnd,dum,dum,dum,dum,cofx,cofy,grhs,usol,idmn,w, 2 pertrb,ierror) c c compute discretization error c err = 0.0 do 107 j=1,ny y = c+float(j-1)*dly do 106 i=1,nx x = a+float(i-1)*dlx err = amax1(err,abs((usol(i,j)-ue(x,y))/ue(x,y))) 106 continue 107 continue err4 = err iw = int(w(1)) print 1001 , ierror,err2,err4,iw c c 1001 format (1h1,20x,25hsubroutine sepeli example/// 1 20x,46hthe output from the ncar control data 7600 was// 2 20x,10hierror = 0/ 3 20x,47hsecond order discretization error = 9.78910e-05/ 4 20x,47hfourth order discretization error = 1.47351e-06/ 5 20x,33hrequired length of w array = 1118// 6 20x,32hthe output from your computer is// 7 20x,8hierror =i2/ 8 20x,36hsecond order discretization error = e12.5/ 9 20x,36hfourth order discretization error = e12.5/ + 20x,29hrequired length of w array = i4) c end