c fishpk41 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 this program illustrates the use of the subroutine pois3d to solve c the equation c c (d/dx)(du/dx) + (d/dy)(du/dy) + (1+z)**2*(d/dz)(du/dz) c c - 2*(1+z)*(du/dz) = 2*sin(x)*sin(y)*(1+z)**4 (1) c c on the parallelepiped -pi .lt. x .lt. pi, -pi .lt. y .lt. pi, c 0 .lt. z .lt. 1 with boundary conditions c c u periodic in x c c u periodic in y c c (du/dz)(x,y,0) = 4*sin(x)*sin(y) -pi .lt. x,y .lt. pi (2) c c u(x,y,1) = 16*sin(x)*sin(y) -pi .lt. x,y .lt. pi (3) c c using a finite difference grid with deltax (= dx) = 2*pi/30 , c deltay (= dy) = 2*pi/30, and deltaz (= dz) = 1/10. c to set up the finite difference equations we define the grid c points c c x(i) = -pi + (i-1)*dx i=1,2,...,31 c c y(j) = -pi + (j-1)*dy j=1,2,...,31 c c z(k) = (k-1)*dz k=1,2,...,11 c c and let v(i,j,k) be an approximation to u(x(i),y(j),z(k)). c numbering the grid points in this fashion gives the set of c unknowns as v(i,j,k) for i=1,2,...,30, j=1,2,...,30, k=1,2,...,10. c hence, in the program l=30, m = 30, and n = 10. at the interior c grid point (x(i),y(j),z(k)), we replace all derivatives in c equation (1) by second order central finite differences and c collect coefficients of v(i,j,k) to get the finite difference c equation c c (v(i-1,j,k) - 2v(i,j,k) + v(i+1,j,k))/dx**2 c c + (v(i,j-1,k) - 2v(i,j,k) + v(i,j+1,k))/dy**2 c c + a(k)v(i,j,k-1) + b(k)v(i,j,k) + c(k)v(i,j,k+1) = f(i,j,k) (4) c c where for k=2,3,...,9 c c a(k) = (1+z(k))**2/dz**2 + (1+z(k))/dz c c b(k) = -2(1+z(k))**2/dz**2 c c c(k) = (1+z(k))**2/dz**2 - (1+z(k))/dz c c f(i,j,k) = 2sin(x(i))*sin(y(j))*(1+z(k))**4 for i,j=1,2,...,30. c c to obtain equations for k=1, we replace the derivative in c equation (2) by a second order central finite difference approx- c imation, use this equation to eliminate the virtual unknown c v(i,j,0) in equation (4) and arrive at the equation c c (v(i-1,j,1) -2v(i,j,1) + v(i+1,j,1))/dx**2 c c + (v(i,j-1,1) -2v(i,j,1) + v(i,j+1,1))/dy**2 c c + b(1)v(i,j,1) + c(1)v(i,j,2) = f(i,j,1) c c where c b(1) = -c(1) = -2(1+z(1))**2/dz**2 = -2/dz**2 c c f(i,j,1) = (10 + 8/dz)sin(x(i))*sin(y(j)) c c for i,j=1,2,...,30. for completeness we set a(1) = 0. c c to obtain equations for k=10, we incorporate equation (3) into c equation (4) by setting c c v(i,j,11) = u(x(i),y(j),1) = 16sin(x(i))*sin(y(j)) c c and arrive at the equation c c (v(i-1,j,10) - 2v(i,j,10) + v(i+1,j,10))/dx**2 c c + (v(i,j-1,10) - 2v(i,j,10) + v(i,j+1,10))/dy**2 c c + a(10)v(i,j,9) + b(10)v(i,j,10) = f(i,j,10) c c where c c a(10) = (1+z(10))**2/dz**2 + (1+z(10))/dz c c b(10) = -2(1+z(10))**2/dz**2 c c f(i,j,10) = 2sin(x(i))*sin(y(j))*((1+z(10))**4 c -8*((1+z(10))**2/dz**2 - (1+z(10))/dz)) c c for i,j=1,2,...,30. c c for completeness, we set c(10) = 0. hence, in the program, c nperod = 1. c the periodicity conditions on u give the conditions c c v(0,j,k) = v(30,j,k) and v(31,j,k) = v(1,j,k) c for j=1,2,...,30 and k=1,2,...,10, c and c v(i,0,k) = v(i,30,k) and v(i,31,k) = v(i,1,k) c for i=1,2,...,30 and k=1,2,...,10. c c hence, in the program lperod = mperod = 0. c dimension f(32,33,10),a(10) ,b(10) ,c(10) , 1 w(350) ,x(30) ,y(30) ,z(10) c c from the dimension statement we get that ldimf = 32, mdimf = 33, c and note that w has been dimensioned according to its description. c ldimf = 32 mdimf = 33 pi = pimach(dum) lperod = 0 l = 30 dx = 2.*pi/float(l) c1 = 1./dx**2 mperod = 0 m = 30 dy = 2.*pi/float(m) c2 = 1./dy**2 nperod = 1 n = 10 dz = 1./float(n) dzsq = 1./dz**2 c c generate grid points for later use. c do 101 i=1,l x(i) = -pi+float(i-1)*dx 101 continue do 102 j=1,m y(j) = -pi+float(j-1)*dy 102 continue c c generate coefficients c a(1) = 0. b(1) = -2.*dzsq c(1) = -b(1) z(1) = 0. do 103 k=2,n z(k) = float(k-1)*dz t = 1.+z(k) a(k) = t**2*dzsq+t/dz b(k) = -2.*t**2*dzsq c(k) = t**2*dzsq-t/dz 103 continue c c generate right side of equation c do 106 i=1,l do 105 j=1,m do 104 k=2,n f(i,j,k) = 2.*sin(x(i))*sin(y(j))*(1.+z(k))**4 104 continue 105 continue 106 continue do 108 i=1,l do 107 j=1,l f(i,j,1) = (10.+8./dz)*sin(x(i))*sin(y(j)) f(i,j,n) = f(i,j,n)-c(n)*16.*sin(x(i))*sin(y(j)) 107 continue 108 continue c(n) = 0. c c call pois3d to solve equations. c call pois3d (lperod,l,c1,mperod,m,c2,nperod,n,a,b,c,ldimf,mdimf, 1 f,ierror,w) c c compute discretization error. the exact solution is c c u(x,y,z) = sin(x)*sin(y)*(1+z)**4 c err = 0. do 111 i=1,l do 110 j=1,m do 109 k=1,n t = abs(f(i,j,k)-sin(x(i))*sin(y(j))*(1.+z(k))**4) if (t .gt. err) err = t 109 continue 110 continue 111 continue print 1001 , ierror,err stop c 1001 format (1h1,20x,25hsubroutine pois3d example/// 1 10x,46hthe output from the ncar control data 7600 was// 2 32x,10hierror = 0/ 3 18x,34hdiscretization error = 2.93277e-02// 4 10x,32hthe output from your computer is// 5 32x,8hierror =,i2/18x,22hdiscretization error =,e12.5) c end c compute discretization error. the exact solution is