! D03RBF Example Program Text ! Mark 24 Release. NAG Copyright 2012. Module d03rbfe_mod ! D03RBF Example Program Module: ! Parameters and User-defined Routines ! .. Use Statements .. Use nag_library, Only: nag_wp ! .. Implicit None Statement .. Implicit None ! .. Parameters .. Real (Kind=nag_wp), Parameter :: one = 1.0_nag_wp Real (Kind=nag_wp), Parameter :: zero = 0.0_nag_wp Real (Kind=nag_wp), Parameter :: twant(2) = (/0.25_nag_wp,one/) Integer, Parameter :: itrace = 0, nin = 5, nout = 6, & npde = 2 ! .. Local Scalars .. Integer :: iout Contains Subroutine pdeiv(npts,npde,t,x,y,u) ! .. Parameters .. Real (Kind=nag_wp), Parameter :: eps = 0.001_nag_wp ! .. Scalar Arguments .. Real (Kind=nag_wp), Intent (In) :: t Integer, Intent (In) :: npde, npts ! .. Array Arguments .. Real (Kind=nag_wp), Intent (Out) :: u(npts,npde) Real (Kind=nag_wp), Intent (In) :: x(npts), y(npts) ! .. Local Scalars .. Real (Kind=nag_wp) :: a Integer :: i ! .. Intrinsic Procedures .. Intrinsic :: exp ! .. Executable Statements .. Do i = 1, npts a = (4.0_nag_wp*(y(i)-x(i))-t)/(32.0_nag_wp*eps) If (a<=zero) Then u(i,1) = 0.75_nag_wp - 0.25_nag_wp/(one+exp(a)) u(i,2) = 0.75_nag_wp + 0.25_nag_wp/(one+exp(a)) Else a = -a u(i,1) = 0.75_nag_wp - 0.25_nag_wp*exp(a)/(one+exp(a)) u(i,2) = 0.75_nag_wp + 0.25_nag_wp*exp(a)/(one+exp(a)) End If End Do Return End Subroutine pdeiv Subroutine inidom(maxpts,xmin,xmax,ymin,ymax,nx,ny,npts,nrows,nbnds, & nbpts,lrow,irow,icol,llbnd,ilbnd,lbnd,ierr) ! .. Use Statements .. Use nag_library, Only: d03ryf ! .. Scalar Arguments .. Real (Kind=nag_wp), Intent (Out) :: xmax, xmin, ymax, ymin Integer, Intent (Inout) :: ierr Integer, Intent (In) :: maxpts Integer, Intent (Out) :: nbnds, nbpts, npts, nrows, nx, & ny ! .. Array Arguments .. Integer, Intent (Inout) :: icol(*), ilbnd(*), irow(*), & lbnd(*), llbnd(*), lrow(*) ! .. Local Scalars .. Integer :: i, ifail, j, leniwk ! .. Local Arrays .. Integer :: icold(105), ilbndd(28), & irowd(11), iwk(122), & lbndd(72), llbndd(28), lrowd(11) Character (33) :: pgrid(11) ! .. Executable Statements .. icold(1:105) = (/0,1,2,0,1,2,3,4,5,6,7,8,9,10,0,1,2,3,4,5,6,7,8,9,10, & 0,1,2,3,4,5,6,7,8,9,10,0,1,2,3,4,5,8,9,10,0,1,2,3,4,5,6,7,8,9,10,0, & 1,2,3,4,5,6,7,8,9,10,0,1,2,3,4,5,6,7,8,9,10,0,1,2,3,4,5,6,7,8,0,1,2, & 3,4,5,6,7,8,0,1,2,3,4,5,6,7,8/) ilbndd(1:28) = (/1,2,3,4,1,4,1,2,3,4,3,4,1,2,12,23,34,41,14,41,12,23, & 34,41,43,14,21,32/) irowd(1:11) = (/0,1,2,3,4,5,6,7,8,9,10/) lbndd(1:72) = (/2,4,15,26,37,46,57,68,79,88,98,99,100,101,102,103,104, & 96,86,85,84,83,82,70,59,48,39,28,17,6,8,9,10,11,12,13,18,29,40,49, & 60,72,73,74,75,76,77,67,56,45,36,25,33,32,42,52,53,43,1,97,105,87, & 81,3,7,71,78,14,31,51,54,34/) llbndd(1:28) = (/1,2,11,18,19,24,31,37,42,48,53,55,56,58,59,60,61,62, & 63,64,65,66,67,68,69,70,71,72/) lrowd(1:11) = (/1,4,15,26,37,46,57,68,79,88,97/) nx = 11 ny = 11 ! Check MAXPTS against rough estimate of NPTS npts = nx*ny If (maxpts1. If (iout/=2 .Or. level/=1) Cycle levels ! Get exact solution Call pdeiv(npts,npde,t,x,y,uex) Write (nout,*) Write (nout,99999) t Write (nout,*) Write (nout,99998) 'x', 'y', 'approx u', 'exact u', 'approx v', & 'exact v' Write (nout,*) ipsol = lsol(level) Do i = 1, npts, 2 Write (nout,99997) x(i), y(i), sol(ipsol+i), uex(i,1), & sol(ipsol+npts+i), uex(i,2) End Do Write (nout,*) End Do levels Return 99999 Format (' Solution at every 2nd grid point in level 1 at time ',F8.4, & ':') 99998 Format (7X,A,10X,A,8X,A,5X,A,4X,A,4X,A) 99997 Format (6(1X,E11.4)) End Subroutine monitr End Module d03rbfe_mod Program d03rbfe ! D03RBF Example Main Program ! .. Use Statements .. Use nag_library, Only: d03rbf Use d03rbfe_mod, Only: bndary, inidom, iout, itrace, monitr, nag_wp, & nin, nout, npde, one, pdedef, pdeiv, twant, zero ! .. Implicit None Statement .. Implicit None ! .. Local Scalars .. Real (Kind=nag_wp) :: tols, tolt, tout, ts Integer :: i, ifail, ind, j, leniwk, & lenlwk, lenrwk, maxlev, mxlev, & npts ! .. Local Arrays .. Real (Kind=nag_wp) :: dt(3) Real (Kind=nag_wp), Allocatable :: optr(:,:), rwk(:) Integer, Allocatable :: iwk(:) Integer :: opti(4) Logical, Allocatable :: lwk(:) ! .. Executable Statements .. Write (nout,*) 'D03RBF Example Program Results' ! Skip heading in data file Read (nin,*) Read (nin,*) npts, mxlev leniwk = npts*(5*mxlev+14) + 2 + 7*mxlev lenlwk = 2*npts lenrwk = npts*npde*(5*mxlev+9+18*npde) + 2*npts Allocate (rwk(lenrwk),iwk(leniwk),lwk(lenlwk),optr(3,npde)) ind = 10 ts = zero dt(1) = 0.001_nag_wp dt(2) = 1.0E-7_nag_wp dt(3) = zero tols = 0.1_nag_wp tolt = 0.05_nag_wp opti(1) = 5 maxlev = opti(1) opti(2:4) = 0 optr(1:3,1:npde) = one ! Call main routine Do iout = 1, 2 tout = twant(iout) ! ifail: behaviour on error exit ! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft ifail = 0 Call d03rbf(npde,ts,tout,dt,tols,tolt,inidom,pdedef,bndary,pdeiv, & monitr,opti,optr,rwk,lenrwk,iwk,leniwk,lwk,lenlwk,itrace,ind,ifail) ! Print statistics Write (nout,99999) 'Statistics:' Write (nout,99998) 'Time = ', ts Write (nout,99997) 'Total number of accepted timesteps =', iwk(1) Write (nout,99997) 'Total number of rejected timesteps =', iwk(2) Write (nout,99996) maxlev = opti(1) Do j = 1, maxlev If (iwk(j+2)/=0) Then Write (nout,'(I6,4I10)') j, (iwk(j+2+i*maxlev),i=0,3) End If End Do Write (nout,99995) Do j = 1, maxlev If (iwk(j+2)/=0) Then Write (nout,'(I6,2I14)') j, (iwk(j+2+i*maxlev),i=4,5) End If End Do Write (nout,*) End Do 99999 Format (1X,A) 99998 Format (1X,A,F8.4) 99997 Format (1X,A,I5) 99996 Format (1X/' T o t a l n u m b e r o f '/' ', & ' Residual Jacobian Newton Lin sys'/' ', & ' evals evals iters iters'/' At level ') 99995 Format (1X/' M a x i m u m n u m b e r o f'/' ', & ' Newton iters Lin sys iters '/' At level ') End Program d03rbfe