! fishprb.f 24 February 1998 ! program fishprb ! !*********************************************************************** ! !! FISHPRB calls the FISHPACK tests. ! ! KPRINT = 0 Quick checks - No printing. ! Driver - Short pass or fail message printed. ! 1 Quick checks - No message printed for passed tests, ! short message printed for failed tests. ! Driver - Short pass or fail message printed. ! 2 Quick checks - Print short message for passed tests, ! fuller information for failed tests. ! Driver - Pass or fail message printed. ! 3 Quick checks - Print complete quick check results. ! Driver - Pass or fail message printed. ! integer ipass integer kprint integer nfail ! nfail = 0 kprint=2 call xermax(1000) if (kprint .le. 1) then call xsetf(0) else call xsetf(1) end if ! ! Test HWSCRT ! call test01(kprint,ipass) if (ipass .eq. 0) nfail = nfail + 1 ! ! Test HWSPLR ! call test02(kprint,ipass) if (ipass .eq. 0) nfail = nfail + 1 ! ! Test HWSCYL ! call test03(kprint,ipass) if (ipass .eq. 0) nfail = nfail + 1 ! ! Test HWSSSP ! call test04(kprint,ipass) if (ipass .eq. 0) nfail = nfail + 1 ! ! Test HWSCSP ! call test05(kprint,ipass) if (ipass .eq. 0) nfail = nfail + 1 call test06(kprint,ipass) if ( ipass .eq. 0 ) nfail = nfail + 1 ! ! Test GENBUN ! call test07(kprint,ipass) if (ipass .eq. 0) nfail = nfail + 1 ! ! Test BLKTRI ! call test08(kprint,ipass) if (ipass .eq. 0) then nfail = nfail + 1 end if ! ! Test HWSCRT ! call test09 ! ! Test HWSPLR ! call test10 ! ! Test BLKTRI ! call test11 ! ! Write PASS or FAIL message ! write(*,*)' ' if (nfail .eq. 0) then write(*,*)'FISHPRB passed all tests' else write(*,*)'FISHPRB - Warning!' write(*,*)' Failed ',nfail,' tests.' end if stop end subroutine test01(kprint, ipass) ! !*********************************************************************** ! !! TEST01 illustrates the use of HWSCRT. ! ! The equation solved is: ! ! (D/DX)(DU/DX) + (D/DY)(DU/DY) - 4*U ! ! = (2 - (4 + PI**2/4)*X**2)*COS((Y+1)*PI/2) ! ! on the rectangle 0 .LT. X .LT. 2, -1 .LT. Y .LT. 3 with the ! with the boundary conditions ! ! U(0,Y) = 0 ! -1 .LE. Y .LE. 3 ! (DU/DX)(2,Y) = 4*COS((Y+1)*PI/2) ! ! and with U periodic in Y. ! ! The X-interval will be divided into 40 intervals and the ! Y-interval will be divided into 80 intervals. ! real a real b real bda real bdb(81) real bdc real bdd real c real d real elmbda real ermax real err real f(45,82) integer i integer idimf integer ierror integer ipass integer j integer kprint integer m integer mbdcnd integer n integer nbdcnd real pertrb real pi real piby2 real pimach real pisq real w(1200) real x(41) real y(81) real z ! ! From dimension statement we get value of IDIMF. Also note that W ! is dimensioned 6*(N+1) + 8*(M+1). ! idimf = 45 ermax=0.001 a = 0.0 b = 2.0 m = 40 mbdcnd = 2 c = -1.0 d = 3.0 n = 80 nbdcnd = 0 elmbda = -4.0 ! ! Auxiliary quantities. ! pi = pimach() piby2 = pi/2. pisq = pi**2 ! ! Generate and store grid points for the purpose of computing ! boundary data and the right side of the Helmholtz equation. ! do i=1,m+1 x(i) = (i-1)/20.0 end do do j=1,n+1 y(j) = -1.0+(j-1)/20.0 end do ! ! Generate boundary data. ! do j=1,n+1 bdb(j) = 4.0*cos((y(j)+1.)*piby2) end do ! ! BDA, BDC, and BDD are dummy variables. ! do j=1,n+1 f(1,j) = 0.0 end do ! ! Generate right side of equation. ! do i=2,m+1 do j=1,n+1 f(i,j) = (2.-(4.+pisq/4.)*x(i)**2)*cos((y(j)+1.)*piby2) end do end do call hwscrt(a,b,m,mbdcnd,bda,bdb,c,d,n,nbdcnd,bdc,bdd,elmbda,f, & idimf,pertrb,ierror,w) ! ! Compute discretization error. The exact solution is ! U(X,Y) = X**2*COS((Y+1)*PIBY2) ! err = 0.0 do i=1,m+1 do j=1,n+1 z = abs(f(i,j)-x(i)**2*cos((y(j)+1.)*piby2)) if (z .gt. err) err = z end do end do ipass = 1 if (err.gt.ermax) ipass = 0 if (kprint.eq.0) return if (kprint.ge.2 .or. ipass.eq.0) then write(*,*)' ' write(*,1001) ierror,err,int(w(1)) if (ipass.eq.1) then write(*,*)'Test passed.' else write(*, *)'Test failed.' end if end if return 1001 format (' TEST01: HWSCRT example'/ & ' the output from the NCAR Control Data 7600 was'/, & ' ierror = 0'/, & ' discretization error = 5.36508e-04'/, & ' required length of w array = 880'/, & ' the output from your computer is'/, & ' ierror =',i2,/, & ' discretization error =',1pe12.5,/, & ' required length of w array =',i4) end subroutine test02(kprint, ipass) ! !*********************************************************************** ! !! TEST02 illustrates the use of HWSPLR. ! ! The equation solved is ! ! (1/R)(D/DR)(R*(DU/DR)) + (1/R**2)(D/DTHETA)(DU/DTHETA) = 16*R**2 ! ! on the quarter-disk 0 .LT. R .LT. 1, 0 .LT. THETA .LT. PI/2 ! with the boundary conditions ! ! U(1,THETA) = 1 - COS(4*THETA), 0 .LE. THETA .LE. 1 ! ! and ! ! (DU/DTHETA)(R,0) = (DU/DTHETA)(R,PI/2) = 0, 0 .LE. R .LE. 1. ! ! (Note that the solution U is unspecified at R = 0.) ! The R-interval will be divided into 50 panels and the ! THETA-interval will be divided into 48 panels. ! real a real b real bda real bdb real bdc(51) real bdd(51) real c real d real elmbda real ermax real err real f(100,50) integer i integer idimf integer ierror integer ipass integer j integer kprint integer m integer mbdcnd integer n integer nbdcnd real pertrb real pi real pimach real r(51) real theta(49) real w(1200) real z ! ! From dimension statement we get value of IDIMF. Also note that W ! is dimensioned 6*(N+1) + 8*(M+1). ! idimf = 100 ermax=0.001 a = 0. b = 1. m = 50 mbdcnd = 5 c = 0. pi = pimach() d = pi/2.0 n = 48 nbdcnd = 3 elmbda = 0. ! ! Generate and store grid points for the purpose of computing ! boundary data and the right side of the Poisson equation. ! do i=1,m+1 r(i) = (i-1)/50.0 end do do j=1,n+1 theta(j) = (j-1)*pi/96.0 end do ! ! Generate boundary data. ! do i=1,m+1 bdc(i) = 0.0 bdd(i) = 0.0 end do ! ! BDA and BDB are dummy variables. ! do j=1,n+1 f(m+1,j) = 1.0-cos(4.0*theta(j)) end do ! ! Generate right side of equation. ! do i=1,m do j=1,n+1 f(i,j) = 16.*r(i)**2 end do end do call hwsplr(a,b,m,mbdcnd,bda,bdb,c,d,n,nbdcnd,bdc,bdd,elmbda,f, & idimf,pertrb,ierror,w) ! ! Compute discretization error. The exact solution is ! U(R,THETA) = R**4*(1 - COS(4*THETA)) ! err = 0. do i=1,m+1 do j=1,n+1 z = abs(f(i,j)-r(i)**4*(1.-cos(4.*theta(j)))) if (z .gt. err) err = z end do end do ipass = 1 if (err.gt.ermax) ipass = 0 if (kprint.eq.0) return if (kprint.ge.2 .or. ipass.eq.0) then write(*,*)' ' write(*,1001) ierror,err,int(w(1)) if (ipass.eq.1) then write(*,*)'Test passed.' else write(*, *)'Test failed.' end if end if return 1001 format (' TEST02: HWSPLR example'/, & ' the output from the NCAR Control Data 7600 was',/, & ' ierror = 0',/, & ' discretization error = 6.19134e-04',/, & ' required length of w array = 882',/, & ' the output from your computer is',/, & ' ierror =',i2,/, & ' discretization error =',1pe12.5,/, & ' required length of W array =',i4) end subroutine test03(kprint, ipass) ! !*********************************************************************** ! !! TEST03 illustrates the use of HWSCYL. ! ! The equation solved is ! ! (1/R)(D/DR)(R*(DU/DR)) + (D/DZ)(DU/DZ) ! ! = (2*R*Z)**2*(4*Z**2 + 3*R**2) ! ! on the rectangle 0 .LT. R .LT. 1, 0 .LT. Z .LT. 1 with the ! boundary conditions ! ! U(0,Z) unspecified ! 0 .LE. Z .LE. 1 ! (DU/DR)(1,Z) = 4*Z**4 ! ! and ! ! (DU/DZ)(R,0) = 0 ! 0 .LE. R .LE. 1 ! (DU/DZ)(R,1) = 4*R**4 . ! ! The R-interval will be divided into 50 panels and the ! Z-interval will be divided into 100 panels. ! real a real b real bda(101) real bdb(101) real bdc(51) real bdd(51) real c real d real elmbda real ermax real err real f(75,105) integer i integer idimf integer ierror integer ipass integer j integer kprint integer m integer mbdcnd integer n integer nbdcnd real pertrb real r(51) real w(1200) real x real z(101) ! ! From dimension statement we get value of IDIMF. Also note that W ! is dimensioned 6*(N+1) + 8*(M+1). ! idimf = 75 ermax=0.001 a = 0. b = 1. m = 50 mbdcnd = 6 c = 0. d = 1. n = 100 nbdcnd = 3 elmbda = 0. ! ! Generate and store grid points for the purpose of computing ! boundary data and the right side of the Poisson equation. ! do i=1,m+1 r(i) = (i-1)/50.0 end do do j=1,n+1 z(j) = (j-1)/100.0 end do ! ! Generate boundary data. ! do j=1,n+1 bdb(j) = 4.*z(j)**4 end do do i=1,m+1 bdc(i) = 0. bdd(i) = 4.*r(i)**4 end do ! ! BDA is a dummy variable. ! ! ! Generate right side of equation. ! do i=1,m+1 do j=1,n+1 f(i,j) = 4.*r(i)**2*z(j)**2*(4.*z(j)**2+3.*r(i)**2) end do end do call hwscyl(a,b,m,mbdcnd,bda,bdb,c,d,n,nbdcnd,bdc,bdd,elmbda,f, & idimf,pertrb,ierror,w) ! ! Compute discretization error by minimizing over all A the function ! NORM(F(I,J) - A*1 - U(R(I),Z(J))). The exact solution is ! U(R,Z) = (R*Z)**4 + arbitrary constant. ! x = 0.0 do i=1,m+1 do j=1,n+1 x = x+f(i,j)-(r(i)*z(j))**4 end do end do x = x/((n+1)*(m+1)) do i=1,m+1 do j=1,n+1 f(i,j) = f(i,j)-x end do end do err = 0.0 do i=1,m+1 do j=1,n+1 x = abs(f(i,j)-(r(i)*z(j))**4) if (x .gt. err) err = x end do end do ipass = 1 if (err.gt.ermax) ipass = 0 if (kprint.eq.0) return if (kprint.ge.2 .or. ipass.eq.0) then write(*,*)' ' write(*,1001) ierror,pertrb,err,int(w(1)) if (ipass.eq.1) then write(*,*)'Test passed.' else write(*, *)'Test failed.' end if end if return 1001 format (' TEST03: HWSCYL example'/ & ' the output from the NCAR Control Data 7600 was',/, & ' ierror = 0',/, & ' pertrb = 2.26734e-04',/, & ' discretization error = 3.73672e-04',/, & ' required length of w array = 1118',/, & ' the output from your computer is',/, & ' ierror =',i2,/, & ' pertrb =',e12.5,/, & ' discretization error =',1pe12.5,/, & ' required length of w array =',i4) end subroutine test04(kprint, ipass) ! !*********************************************************************** ! !! TEST04 illustrates the use of HWSSSP. ! real bdpf real bdps real bdtf(73) real bdts real dphi real dtheta real elmbda real ermax real err real f(19,73) integer i integer idimf integer ierror integer ipass integer j integer kprint integer m integer mbdcnd integer n integer nbdcnd real pertrb real pf real pi real pimach real ps real sinp(73) real sint(19) real tf real ts real w(1200) real z ! ! The value of IDIMF is the first dimension of F. W is ! dimensioned 11*(M+1)+6*(N+1)=647 since M=18 and N=72. ! pi = pimach() ermax=0.005 ts = 0.0 tf = pi/2. m = 18 mbdcnd = 6 ps = 0.0 pf = pi+pi n = 72 nbdcnd = 0 elmbda = 0. idimf = 19 ! ! Generate sines for use in subsequent computations ! dtheta = tf/m do i=1,m+1 sint(i) = sin((i-1)*dtheta) end do dphi = (pi+pi)/n do j=1,n+1 sinp(j) = sin((j-1)*dphi) end do ! ! Compute right side of equation and store in F ! do j=1,n+1 do i=1,m+1 f(i,j) = 2.-6.*(sint(i)*sinp(j))**2 end do end do ! ! Store derivative data at the equator ! do j=1,n+1 bdtf(j) = 0. end do call hwsssp(ts,tf,m,mbdcnd,bdts,bdtf,ps,pf,n,nbdcnd,bdps,bdpf, & elmbda,f,idimf,pertrb,ierror,w) ! ! Compute discretization error. Since problem is singular, the ! solution must be normalized. ! err = 0.0 do j=1,n+1 do i=1,m+1 z = abs(f(i,j)-(sint(i)*sinp(j))**2-f(1,1)) if (z .gt. err) err = z end do end do ipass = 1 if (err.gt.ermax) ipass = 0 if (kprint.eq.0) return if (kprint.ge.2 .or. ipass.eq.0) then write(*,*)' ' write(*,1001) ierror,err,int(w(1)) if (ipass.eq.1) then write(*,*)'Test passed.' else write(*, *)'Test failed.' end if end if return 1001 format (' TEST04: HWSSP example'/, & ' the output from the NCAR Control Data 7600 was',/, & ' ierror = 0',/, & ' discretization error = 3.38107e-03',/, & ' required length of w array = 600',/, & ' the output from your computer is',/, & ' ierror =',i2,/, & ' discretization error =',1pe12.5, /, & ' required length of W array =',i4) end subroutine test05(kprint, ipass) ! !*********************************************************************** ! !! TEST05 illustrates the use of HWSCSP. ! real bdrf real bdrs real bdtf(33) real bdts real ci4 real dr real dtheta real elmbda real ermax real err real f(48,33) integer i integer idimf integer ierror integer intl integer ipass integer j integer kprint integer m integer mbdcnd integer n integer nbdcnd real pertrb real pi real pimach real r(33) real rf real rs real tf real theta(48) real ts real w(1200) real z ! ! The value of IDIMF is the first dimension of F. Since M=36, N=32, ! L=N therefore K=5 and W is dimensioned 2*(L+1)*(K-1) + 6*(M+N) ! + MAX(4*N,6*M) + 14 = 902. ! ermax=0.001 pi = pimach() intl = 0 ts = 0. tf = pi/2. m = 36 mbdcnd = 6 rs = 0. rf = 1. n = 32 nbdcnd = 5 elmbda = 0. idimf = 48 ! ! Generate and store grid points for the purpose of computing the ! boundary data and the right side of the equation. ! dtheta = tf/m do i=1,m+1 theta(i) = (i-1)*dtheta end do dr = 1.0/n do j=1,n+1 r(j) = (j-1)*dr end do ! ! Generate normal derivative data at equator ! do j=1,n+1 bdtf(j) = 0. end do ! ! Compute boundary data on the surface of the sphere ! do i=1,m+1 f(i,n+1) = cos(theta(i))**4 end do ! ! Compute right side of equation ! do i=1,m+1 ci4 = 12.0*cos(theta(i))**2 do j=1,n f(i,j) = ci4*r(j)**2 end do end do call hwscsp(intl,ts,tf,m,mbdcnd,bdts,bdtf,rs,rf,n,nbdcnd,bdrs, & bdrf,elmbda,f,idimf,pertrb,ierror,w) ! ! Compute discretization error ! err = 0. do i=1,m+1 ci4 = cos(theta(i))**4 do j=1,n z = abs(f(i,j)-ci4*r(j)**4) if (z .gt. err) err = z end do end do ipass = 1 if (err.gt.ermax) ipass = 0 if (kprint.ne.0) then if (kprint.ge.2 .or. ipass.eq.0) then write(*,*)' ' write(*,1001) ierror,err,int(w(1)) if (ipass.eq.1) then write(*, *)'Test passed.' else write(*, *)'Test failed.' end if end if end if 1001 format (' TEST05: HWSCSP example 1'/ & ' the output from the NCAR Control Data 7600 was',/, & ' ierror = 0',/, & ' discretization error = 7.99842e-04',/, & ' required length of w array = 775',/, & ' the output from your computer is',/, & ' ierror =',i2,/, & ' discretization error =',1pe12.5,/, & ' required length of w array =',i4) return end subroutine test06(kprint, ipass) ! !*********************************************************************** ! !! TEST06 illustrates the use of HWSCSP. ! ! The following program illustrates the use of HWSCSP to solve ! a three dimensional problem which has longitudinal dependence ! real bdrf real bdrs real bdtf(33) real bdts real dphi real dr real dtheta real elmbda real ermax real err real f(48,33) integer i integer idimf integer ierror integer intl integer ipass integer j integer kprint integer m integer mbdcnd integer n integer nbdcnd real pertrb real pi real pimach real r(33) real rf real rs real si real tf real theta(48) real ts real w(1200) real z ! ! The value of IDIMF is the first dimension of F. Since M=36, N=32, ! L=N therefore K=5 and W is dimensioned 2*(L+1)*(K-1) + 6*(M+N) ! + MAX(4*N,6*M) + 14 = 902. ! ermax=0.001 pi = pimach() intl = 0 ts = 0. tf = pi/2. m = 36 mbdcnd = 2 rs = 0. rf = 1. n = 32 nbdcnd = 1 idimf = 48 dphi = pi/72. elmbda = -2.0*(1.0-cos(dphi))/dphi**2 ! ! Generate and store grid points for the purpose of computing the ! boundary data and the right side of the equation. ! dtheta = tf/m do i=1,m+1 theta(i) = (i-1)*dtheta end do dr = 1.0/n do j=1,n+1 r(j) = (j-1)*dr end do ! ! Generate normal derivative data at equator ! do j=1,n+1 bdtf(j) = 0. end do ! ! Compute boundary data on the surface of the sphere ! do i=1,m+1 f(i,n+1) = sin(theta(i)) end do ! ! Compute right side of the equation ! do j=1,n do i=1,m+1 f(i,j) = 0.0 end do end do call hwscsp(intl,ts,tf,m,mbdcnd,bdts,bdtf,rs,rf,n,nbdcnd,bdrs, & bdrf,elmbda,f,idimf,pertrb,ierror,w) ! ! Compute discretization error (Fourier coefficients) ! err = 0. do i=1,m+1 si = sin(theta(i)) do j=1,n+1 z = abs(f(i,j)-r(j)*si) if (z .gt. err) err = z end do end do if (err.gt.ermax) ipass = 0 if (kprint.eq.0) return if (kprint.ge.2 .or. ipass.eq.0) then write(*,*)' ' write(*,1001) ierror,err,int(w(1)) if (ipass.eq.1) then write(*,*)'Test passed.' else write(*, *)'Test failed.' end if end if return 1001 format (' TEST06: HWSCSP example 2'/ & ' the output from the NCAR Control Data 7600 was',/ & ' ierror = 0',/ & ' discretization error = 5.86824e-05',/ & ' required length of w array = 775',/ & ' the output from your computer is'/, & ' ierror =',i2,/, & ' discretization error =',1pe12.5,/, & ' required length of w array =',i4) end subroutine test07(kprint, ipass) ! !*********************************************************************** ! !! TEST07 illustrates the use of GENBUN. ! real a(20) real b(20) real c(20) real deltax real deltay real dysq real ermax real err real f(25,130) integer i integer idimy integer ierror integer ipass integer j integer kprint integer m integer mperod integer n integer nperod real pi real pimach real s real t real w(1200) real x(20) real y(120) real z ! ! From dimension statement we get value of IDIMY. Also note that ! W(.) is dimensioned 6*N + 5*M. ! ermax=0.01 idimy = 25 mperod = 1 m = 20 deltax = 1.0/m nperod = 0 n = 120 pi = pimach() deltay = 2.0*pi/n ! ! Set the grid points. ! do i=1,m x(i) = real(i-1)/real(m) end do do j=1,n y(j) = -pi + 2*(j-1)*pi/real(n) end do ! ! Generate coefficients. ! s = (deltay/deltax)**2 t = s*deltax a(1) = 0.0 b(1) = -2.0*s c(1) = 2.0*s do i=2,m a(i) = (1.0+x(i))**2*s + (1.0+x(i))*t c(i) = (1.0+x(i))**2*s - (1.0+x(i))*t b(i) = -2.0*(1.0+x(i))**2*s end do c(m) = 0.0 ! ! Generate right side of equation for I = 1 showing introduction of ! boundary data. ! dysq = deltay**2 do j=1,n f(1,j) = dysq*(11. + 8./deltax)*sin(y(j)) end do ! ! Generate the right hand side. ! do i=2,m-1 do j=1,n f(i,j) = dysq*3.0*(1.0+x(i))**4*sin(y(j)) end do end do ! ! Generate right side for I=M showing introduction of ! boundary data. ! do j=1,n f(m,j)=dysq*(3.*(1.+x(m))**4 - 16.*((1.+x(m))/deltax)**2 & +16.*(1.+x(m))/deltax)*sin(y(j)) end do ! ! Call GENBUN to solve the system. ! call genbun(nperod,n,mperod,m,a,b,c,idimy,f,ierror,w) ! ! Compute the discretization error. The exact solution is ! U(X,Y) = (1+X)**4*SIN(Y) ! err = 0.0 do i=1,m do j=1,n z = abs(f(i,j)-(1.+x(i))**4*sin(y(j))) if (z .gt. err) err = z end do end do ipass = 1 if (err.gt.ermax) ipass = 0 if (kprint.eq.0) return if (kprint.ge.2 .or. ipass.eq.0) then write(*,*)' ' write(*,1001) ierror, err, int(w(1)) if (ipass.eq.1) then write(*,*)'Test passed.' else write(*, *)'Test failed.' end if end if return 1001 format (' TEST07: GENBUN example'/ & ' the output from the NCAR Control Data 7600 was',/, & ' ierror = 0',/, & ' discretization error = 7.94113e-03',/, & ' required length of w array = 740',/, & ' the output from your computer is'/, & ' ierror =',i2,/, & ' discretization error =',1pe12.5,/, & ' required length of W array =',i4) end subroutine test08(kprint, ipass) ! !*********************************************************************** ! !! TEST08 illustrates the use of BLKTRI. ! integer m integer n parameter ( m = 50 ) parameter ( n = 63 ) ! integer idimy parameter ( idimy = m ) ! real am(m) real an(n) real bm(m) real bn(n) real cm(m) real cn(n) real deltas real deltat real ermax real err real hds real hdt integer i integer ierror integer iflg integer ipass integer j integer kprint integer mp integer np real s(m) real t(n) real tds real tdt real temp1 real temp2 real temp3 real true real w(1952) real y(idimy,n) ! ermax = 0.001 np = 1 mp = 1 ! ! Generate and store the grid points. ! deltas = 1.0/(m+1) do i = 1, m s(i) = i * deltas end do deltat = 1.0/(n+1) do j = 1, n t(j) = j * deltat end do ! ! Compute the coefficients AM,BM,CM corresponding to the S direction ! hds = deltas/2.0 tds = deltas+deltas do i = 1, m temp1 = 1.0/(s(i)*tds) temp2 = 1.0/((s(i)-hds)*tds) temp3 = 1.0/((s(i)+hds)*tds) am(i) = temp1*temp2 cm(i) = temp1*temp3 bm(i) = -(am(i)+cm(i)) end do ! ! Compute the coefficients AN,BN,CN corresponding to the T direction ! hdt = deltat/2.0 tdt = deltat+deltat do j = 1, n temp1 = 1.0/(t(j)*tdt) temp2 = 1.0/((t(j)-hdt)*tdt) temp3 = 1.0/((t(j)+hdt)*tdt) an(j) = temp1*temp2 cn(j) = temp1*temp3 bn(j) = -(an(j)+cn(j)) end do ! ! Compute the right hand side of the equation. ! do j = 1, n do i = 1, m y(i,j) = 3.75 * s(i) * t(j) * ( s(i)**4 + t(j)**4 ) end do end do ! ! Include nonhomogeneous boundary into right side. ! ! Note that the corner Y(M,N) includes contributions ! from both boundaries. ! do j = 1, n y(m,j) = y(m,j) - cm(m) * t(j)**5 end do do i = 1, m y(i,n) = y(i,n) - cn(n) * s(i)**5 end do ! ! Initialize the code. ! iflg = 0 call blktri(iflg,np,n,an,bn,cn,mp,m,am,bm,cm,idimy,y,ierror,w) if ( ierror .ne. 0 ) then ipass = 0 return end if ! ! Solve the system. ! iflg = 1 call blktri(iflg,np,n,an,bn,cn,mp,m,am,bm,cm,idimy,y,ierror,w) if ( ierror .ne. 0 ) then ipass = 0 return end if ! ! Compute the discretization error. ! err = 0.0 do j = 1, n do i = 1, m true = ( s(i) * t(j) )**5 err = max ( err, abs ( y(i,j) - true ) ) end do end do if ( err .gt. ermax ) then ipass = 0 else ipass = 1 end if ! ! Report the results. ! if ( & ( kprint.eq. 1 .and. ipass.eq.0 ) .or. & kprint .ge. 2 ) then write(*,*)' ' write(*,1001) ierror,err,int(w(1)) if (ipass.eq.1) then write(*, *) 'Test passed' else write(*, *) 'Test failed' end if end if return 1001 format (' TEST08: BLKTRI example'/ & ' the output from the NCAR Control Data 7600 was' /, & ' ierror = 0' /, & ' discretization error = 1.6478e-05' /, & ' required length of w array = 823' /, & ' the output from your computer is' /, & ' ierror =',i2, /, & ' discretization error =',1pe12.5, /, & ' required length of w array =', i4) end subroutine test09 ! !*********************************************************************** ! !! TEST09 tests HWSCRT, 2-D Cartesian coordinate Helmholtz equation: ! ! Uxx + Uyy + Lambda*U = F(X,Y) ! ! We use the simple test case ! ! U = X**2 + Y**2, ! ! for X in (0,1), Y in (0,1). ! integer m integer n integer nwork ! parameter (m=5) parameter (n=5) ! ! NWORK=4*(N+1)+(M+1)*(13+ALOG2(N+1)) ! parameter (nwork=4*(n+1)+(m+1)*(13+3)) ! real a real b real bda(n+1) real bdb(n+1) real bdc(m+1) real bdd(m+1) real c real d real dx real dy real elmbda real f(m+1,n+1) integer i integer idimf integer ierror integer j integer mbdcnd integer nbdcnd real pertrb real w(nwork) real x real y ! write(*,*)' ' write(*,*)'TEST09' write(*,*)'test hwscrt, helmholtz equation,' write(*,*)'2-d cartesian coordinates.' write(*,*)' ' ! ! Define corners of region: ! a=0.0 b=1.0 c=0.0 d=1.0 ! ! Define boundary conditions codes: ! mbdcnd=1 nbdcnd=1 ! ! Set boundary conditions: ! dy=(d-c)/real(n) dx=(b-a)/real(m) x=a do j=1,n+1 y=(j-1)*dy f(1,j)=x*x+y*y end do x=b do j=1,n+1 y=(j-1)*dy f(m+1,j)=x*x+y*y end do y=c j=1 do i=1,m+1 x=(i-1)*dx f(i,j)=x*x+y*y end do y=d do i=1,m+1 x=(i-1)*dx f(i,n+1)=x*x+y*y end do ! ! Set right hand side: ! do i=2,m do j=2,n f(i,j)=4.0 end do end do ! ! Specify lambda. ! elmbda=0.0 idimf=m+1 call hwscrt(a,b,m,mbdcnd,bda,bdb,c,d,n,nbdcnd,bdc,bdd,elmbda,f, & idimf,pertrb,ierror,w) write(*,*)' ' if ( ierror.ne.0)then write(*,*)'hwscrt returned error flag ierror=',ierror else write(*,*)'hwscrt computed approximate solution.' write(*,*)' ' if ( pertrb.ne.0.0)then write(*,*)'required perturbation to be subtracted=',pertrb write(*,*)' ' end if do i=1,m+1 write(*,'(1x,6g12.4)')(f(i,j),j=1,n+11) end do end if return end subroutine test10 ! !*********************************************************************** ! !! TEST10 tests HWSPLR, 2-D polar coordinate system Helmholtz equation: ! ! (1/R) D/DR R DU/DR + (1/R*R) D/DTHETA DU/DTHETA + Lambda*U = F(X,Y) ! ! We use the simple test case U=R*R, for R in (1,2), Y in (0,PI/2). ! integer m integer n integer nwork ! parameter (m=5) parameter (n=5) ! ! NWORK=4*(N+1)+(M+1)*(13+ALOG2(N+1)) ! parameter (nwork=4*(n+1)+(m+1)*(13+3)) ! real a real b real bda(n+1) real bdb(n+1) real bdc(m+1) real bdd(m+1) real c real d real elmbda real f(m+1,n+1) integer i integer idimf integer ierror integer j integer mbdcnd integer nbdcnd real pertrb real r real w(nwork) ! write(*,*)' ' write(*,*)'TEST10' write(*,*)'test hwsplr, helmholtz equation,' write(*,*)'2-d polar coordinates.' write(*,*)' ' ! ! Define corners of region: ! a=1.0 b=2.0 c=0.0 d=3.14159265/2.0 ! ! Define boundary conditions codes: ! mbdcnd=1 nbdcnd=3 ! ! Set boundary conditions: ! r=a do j=1,n+1 f(1,j)=r*r end do r=b do j=1,n+1 f(m+1,j)=r*r end do do i=1,m+1 bdc(i)=0.0 end do do i=1,m+1 bdd(i)=0.0 end do ! ! Set the right hand side: ! do i=2,m do j=1,n+1 f(i,j)=4.0 end do end do ! ! Specify lambda. ! elmbda=0.0 ! ! IDIMF is the first dimension of F as declared here. ! idimf=m+1 call hwsplr(a,b,m,mbdcnd,bda,bdb,c,d,n,nbdcnd,bdc,bdd,elmbda,f, & idimf,pertrb,ierror,w) write(*,*)' ' if ( ierror.ne.0)then write(*,*)'HWSPLR returned error flag ierror=',ierror else write(*,*)'HWSPLR computed approximate solution.' write(*,*)' ' if ( pertrb.ne.0.0 ) then write(*,*)'required perturbation to be subtracted=',pertrb write(*,*)' ' end if do i=1,m+1 write(*,'(1x,6g12.4)')(f(i,j),j=1,n+1) end do end if return end subroutine test11 ! !*********************************************************************** ! !! TEST11 uses BLKTRI to solve a simple PDE. ! ! The domain is the unit square, 0 <= x,y <= 1. ! ! There are M vertical grid lines, and N horizontal grid lines, ! equally spaced in the interior of the square. For simplicity, ! we take M = N. This does not count the extra lines that define the ! boundaries. ! ! At each interior node, we wish to approximate the solution ! of a partial differential equation: ! ! d2u/dx2 + d2u/dy2 = 4. ! ! Along the boundaries, we apply the Dirichlet condition ! u(x,y) = x**2 + y**2. ! (In fact, this formula yields the exact solution to the PDE over ! the whole square). ! ! Our partial differential equation is approximated by ! ! u(i,j+1) ! + u(i-1,j) - 4.0*u(i,j) + u(i+1,j) ! + u(i,j-1) = 4.0 * h**2 ! ! Here "H" is the grid spacing. ! ! This equation is written at every interior grid point. For ! grid points next to the boundary, it is necessary to evaluate ! any references to U on the boundary, and carry those terms to ! the right hand side. ! ! The FISHPACK routine BLKTRI is used to handle this system. ! integer m parameter ( m = 5 ) ! ! We choose to set N = M to simplify the equations. ! integer n parameter ( n = m ) ! integer idimy parameter ( idimy = m ) ! integer nw parameter (nw = 100) ! real am(m) real an(n) real bm(m) real bn(n) real cm(m) real cn(n) real err real errmax integer i integer ierror integer iflg integer imax integer j integer jmax integer mp integer np real u ! ! Replace this dimension by something more reasonable. ! real w(nw) real x(0:m+1) real y(0:n+1) real yy(idimy,n) ! ! We set MP and NP to indicate that the boundary conditions are ! not periodic. ! np = 1 mp = 1 ! ! This is the spacing between nodes. ! delta = 1.0 / real ( n+1 ) ! ! Generate and store the grid points. ! This is just for our convenience. ! write(*,*)' ' write(*,*)'TEST11' write(*,*)' A simple example of the direct use of BLKTRI.' write(*,*)' ' write(*,*)' ' write(*,*)'I, X(I)' write(*,*)' ' do i = 0, m+1 x(i) = real ( i ) / real ( m+1 ) write(*,*)i,x(i) end do write(*,*)' ' write(*,*)'J, Y(J)' write(*,*)' ' do j = 0, n+1 y(j) = real ( j ) / real ( n+1 ) write(*,*)j,y(j) end do ! ! Compute the coefficients AM, BM, CM corresponding to the X direction. ! do i = 1, m am(i) = 1.0 bm(i) = -2.0 cm(i) = 1.0 end do ! ! Compute the coefficients AN, BN, CN corresponding to the Y direction. ! do j = 1, n an(j) = 1.0 bn(j) = -2.0 cn(j) = 1.0 end do ! ! Set the right hand side of the equations. ! do j = 1, n do i = 1, m yy(i,j) = 4.0 * delta**2 end do end do ! ! Set the boundary conditions implicitly, by carrying terms involving ! boundary values to the right hand side, and resetting the corresponding ! coefficients to 0. ! i = 1 do j = 1, n uleft = x(i-1)**2 + y(j)**2 yy(i,j) = yy(i,j)-am(1)*uleft end do am(1) = 0.0 i = m do j = 1, n urite = x(i+1)**2 + y(j)**2 yy(i,j) = yy(i,j)-cm(m)*urite end do cm(m) = 0.0 j = 1 do i = 1, m ubot = x(i)**2 + y(j-1)**2 yy(i,j) = yy(i,j)-an(1)*ubot end do an(1) = 0.0 j = n do i = 1, m utop = x(i)**2 + y(j+1)**2 yy(i,j) = yy(i,j)-cn(n)*utop end do cn(n)=0.0 write(*,*)' ' write(*,*)'I,J,YY(I,J)' write(*,*)' ' do i = 1, m do j = 1, n write(*,*)i,j,yy(i,j) end do write (*,*)' ' end do ! ! Initialize the code. ! iflg = 0 call blktri ( iflg, np, n, an, bn, cn, mp, m, am, bm, cm, idimy, & yy, ierror, w ) if ( ierror .ne. 0 ) then write ( *, * ) 'Initial BLKTRI returned IERROR = ',ierror return end if write(*,*)' ' write(*,*)'Array W requires at least ',w(1),' entries.' if ( nw .lt. int(w(1)) ) then write(*,*)' ' write(*,*)'User allocated only ',nw,' entries.' stop end if ! ! Solve the system. ! iflg = 1 call blktri ( iflg, np, n, an, bn, cn, mp, m, am, bm, cm, idimy, & yy, ierror, w ) if ( ierror .ne. 0 ) then write ( *, * ) 'Solve BLKTRI returned IERROR = ',ierror return end if ! ! Compute the discretization error. ! errmax = 0.0 imax = 0 jmax = 0 write(*,*)' ' write(*,*)'I,J,TRUE,YY(I,J)' write(*,*)' ' do j = 1, n do i = 1, m u = x(i)**2 + y(j)**2 write(*,*)i,j,u,yy(i,j) err = abs ( yy(i,j) - u ) if ( err .ge. errmax ) then imax = i jmax = j errmax = err end if end do write(*,*)' ' end do write (*,*) ' ' write (*,*) ' Maximum error is ', errmax write (*,*) ' I,J = ',imax, jmax return end