! D03RAF Example Program Text ! Mark 23 Release. NAG Copyright 2011. MODULE d03rafe_mod ! D03RAF 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 :: activ_energy = 20.0_nag_wp REAL (KIND=nag_wp), PARAMETER :: alpha = 50.0_nag_wp REAL (KIND=nag_wp), PARAMETER :: beta = 300.0_nag_wp REAL (KIND=nag_wp), PARAMETER :: diffusion = 0.1_nag_wp REAL (KIND=nag_wp), PARAMETER :: one = 1.0_nag_wp REAL (KIND=nag_wp), PARAMETER :: reaction_rate = 5.0_nag_wp REAL (KIND=nag_wp), PARAMETER :: xmax = one REAL (KIND=nag_wp), PARAMETER :: ymax = one REAL (KIND=nag_wp), PARAMETER :: zero = 0.0_nag_wp REAL (KIND=nag_wp), PARAMETER :: heat_release = one REAL (KIND=nag_wp), PARAMETER :: xmin = zero REAL (KIND=nag_wp), PARAMETER :: ymin = zero INTEGER, PARAMETER :: itrace = 0, nin = 5, nout = 6, & npde1 = 1, npde2 = 2 ! .. Local Scalars .. REAL (KIND=nag_wp), SAVE :: damkohler INTEGER, SAVE :: iout CONTAINS SUBROUTINE pdeiv1(npts,npde,t,x,y,u) ! .. Implicit None Statement .. IMPLICIT NONE ! .. 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) ! .. Executable Statements .. u(1:npts,1:npde) = one RETURN END SUBROUTINE pdeiv1 SUBROUTINE pdedef1(npts,npde,t,x,y,u,ut,ux,uy,uxx,uxy,uyy,res) ! .. Implicit None Statement .. IMPLICIT NONE ! .. Scalar Arguments .. REAL (KIND=nag_wp), INTENT (IN) :: t INTEGER, INTENT (IN) :: npde, npts ! .. Array Arguments .. REAL (KIND=nag_wp), INTENT (OUT) :: res(npts,npde) REAL (KIND=nag_wp), INTENT (IN) :: u(npts,npde), ut(npts,npde), & ux(npts,npde), & uxx(npts,npde), & uxy(npts,npde), & uy(npts,npde), & uyy(npts,npde), x(npts), & y(npts) ! .. Intrinsic Functions .. INTRINSIC exp ! .. Executable Statements .. res(1:npts,1:npde) = ut(1:npts,1:npde) - diffusion*(uxx(1:npts,1: & npde)+uyy(1:npts,1:npde)) - damkohler*(one+heat_release-u(1:npts, & 1:npde))*exp(-activ_energy/u(1:npts,1:npde)) RETURN END SUBROUTINE pdedef1 SUBROUTINE bndry1(npts,npde,t,x,y,u,ut,ux,uy,nbpts,lbnd,res) ! .. Use Statements .. USE nag_library, ONLY : x02ajf ! .. Implicit None Statement .. IMPLICIT NONE ! .. Scalar Arguments .. REAL (KIND=nag_wp), INTENT (IN) :: t INTEGER, INTENT (IN) :: nbpts, npde, npts ! .. Array Arguments .. REAL (KIND=nag_wp), INTENT (INOUT) :: res(npts,npde) REAL (KIND=nag_wp), INTENT (IN) :: u(npts,npde), ut(npts,npde), & ux(npts,npde), uy(npts,npde), & x(npts), y(npts) INTEGER, INTENT (IN) :: lbnd(nbpts) ! .. Local Scalars .. REAL (KIND=nag_wp) :: tol INTEGER :: i, j ! .. Intrinsic Functions .. INTRINSIC abs ! .. Executable Statements .. tol = 10._nag_wp*x02ajf() DO i = 1, nbpts j = lbnd(i) IF (abs(x(j))<=tol) THEN res(j,1:npde) = ux(j,1:npde) ELSE IF (abs(x(j)-one)<=tol) THEN res(j,1:npde) = u(j,1:npde) - one ELSE IF (abs(y(j))<=tol) THEN res(j,1:npde) = uy(j,1:npde) ELSE IF (abs(y(j)-one)<=tol) THEN res(j,1:npde) = u(j,1:npde) - one END IF END DO RETURN END SUBROUTINE bndry1 SUBROUTINE monit1(npde,t,dt,dtnew,tlast,nlev,ngpts,xpts,ypts,lsol,sol, & ierr) ! .. Implicit None Statement .. IMPLICIT NONE ! .. Scalar Arguments .. REAL (KIND=nag_wp), INTENT (IN) :: dt, dtnew, t INTEGER, INTENT (INOUT) :: ierr INTEGER, INTENT (IN) :: nlev, npde LOGICAL, INTENT (IN) :: tlast ! .. Array Arguments .. REAL (KIND=nag_wp), INTENT (IN) :: sol(*), xpts(*), ypts(*) INTEGER, INTENT (IN) :: lsol(nlev), ngpts(nlev) ! .. Local Scalars .. INTEGER :: i, ipsol, level, npts ! .. Executable Statements .. IF (tlast) THEN ! Print solution IF (iout==2) THEN level = nlev WRITE (nout,99999) level, t WRITE (nout,99998) npts = ngpts(level) ipsol = lsol(level) DO i = 1, npts, 4 WRITE (nout,99997) xpts(i), ypts(i), sol(ipsol+i) END DO WRITE (nout,*) END IF END IF RETURN 99999 FORMAT (1X,'Solution at every 4th grid point in level',I10, & ' at time ',F8.4,':') 99998 FORMAT (1X/7X,'x',10X,'y',8X,'approx u'/) 99997 FORMAT (3(1X,E11.4)) END SUBROUTINE monit1 SUBROUTINE pdeiv2(npts,npde,t,x,y,u) ! .. Use Statements .. USE nag_library, ONLY : x01aaf ! .. Implicit None Statement .. IMPLICIT NONE ! .. 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) :: fourpi ! .. Intrinsic Functions .. INTRINSIC sin ! .. Executable Statements .. fourpi = 4.0_nag_wp*x01aaf(fourpi) u(1:npts,1) = 10.0_nag_wp + (16.0_nag_wp*x(1:npts)*(one-x(1:npts))*y & (1:npts)*(one-y(1:npts)))**2 u(1:npts,2) = -one - alpha*x(1:npts)*y(1:npts) - & beta*sin(fourpi*x(1:npts))*sin(fourpi*y(1:npts)) + & 1.0E4_nag_wp*u(1:npts,1) RETURN END SUBROUTINE pdeiv2 SUBROUTINE pdedef2(npts,npde,t,x,y,u,ut,ux,uy,uxx,uxy,uyy,res) ! .. Use Statements .. USE nag_library, ONLY : x01aaf ! .. Implicit None Statement .. IMPLICIT NONE ! .. Scalar Arguments .. REAL (KIND=nag_wp), INTENT (IN) :: t INTEGER, INTENT (IN) :: npde, npts ! .. Array Arguments .. REAL (KIND=nag_wp), INTENT (OUT) :: res(npts,npde) REAL (KIND=nag_wp), INTENT (IN) :: u(npts,npde), ut(npts,npde), & ux(npts,npde), & uxx(npts,npde), & uxy(npts,npde), & uy(npts,npde), & uyy(npts,npde), x(npts), & y(npts) ! .. Local Scalars .. REAL (KIND=nag_wp) :: fourpi ! .. Local Arrays .. REAL (KIND=nag_wp), ALLOCATABLE :: b1(:) ! .. Intrinsic Functions .. INTRINSIC sin ! .. Executable Statements .. fourpi = 4.0_nag_wp*x01aaf(fourpi) ALLOCATE (b1(npts)) b1(1:npts) = one + alpha*x(1:npts)*y(1:npts) + & beta*sin(fourpi*x(1:npts))*sin(fourpi*y(1:npts)) res(1:npts,1) = ut(1:npts,1) - (uxx(1:npts,1)+uyy(1:npts,1)) - & u(1:npts,1)*(b1(1:npts)-u(1:npts,1)-0.5E-6_nag_wp*u(1:npts,2)) res(1:npts,2) = -0.05_nag_wp*(uxx(1:npts,2)+uyy(1:npts,2)) - & u(1:npts,2)*(-b1(1:npts)+1.0E4_nag_wp*u(1:npts,1)-u(1:npts,2)) DEALLOCATE (b1) RETURN END SUBROUTINE pdedef2 SUBROUTINE bndry2(npts,npde,t,x,y,u,ut,ux,uy,nbpts,lbnd,res) ! .. Use Statements .. USE nag_library, ONLY : x02ajf ! .. Implicit None Statement .. IMPLICIT NONE ! .. Scalar Arguments .. REAL (KIND=nag_wp), INTENT (IN) :: t INTEGER, INTENT (IN) :: nbpts, npde, npts ! .. Array Arguments .. REAL (KIND=nag_wp), INTENT (INOUT) :: res(npts,npde) REAL (KIND=nag_wp), INTENT (IN) :: u(npts,npde), ut(npts,npde), & ux(npts,npde), uy(npts,npde), & x(npts), y(npts) INTEGER, INTENT (IN) :: lbnd(nbpts) ! .. Local Scalars .. REAL (KIND=nag_wp) :: tol INTEGER :: i, j ! .. Intrinsic Functions .. INTRINSIC abs ! .. Executable Statements .. tol = 10.0_nag_wp*x02ajf() DO i = 1, nbpts j = lbnd(i) IF (abs(x(j))<=tol .OR. abs(x(j)-one)<=tol) THEN res(j,1:npde) = ux(j,1:npde) ELSE IF (abs(y(j))<=tol .OR. abs(y(j)-one)<=tol) THEN res(j,1:npde) = uy(j,1:npde) END IF END DO RETURN END SUBROUTINE bndry2 SUBROUTINE monit2(npde,t,dt,dtnew,tlast,nlev,ngpts,xpts,ypts,lsol,sol, & ierr) ! .. Implicit None Statement .. IMPLICIT NONE ! .. Scalar Arguments .. REAL (KIND=nag_wp), INTENT (IN) :: dt, dtnew, t INTEGER, INTENT (INOUT) :: ierr INTEGER, INTENT (IN) :: nlev, npde LOGICAL, INTENT (IN) :: tlast ! .. Array Arguments .. REAL (KIND=nag_wp), INTENT (IN) :: sol(*), xpts(*), ypts(*) INTEGER, INTENT (IN) :: lsol(nlev), ngpts(nlev) ! .. Local Scalars .. INTEGER :: i, ipsol, level, npts ! .. Executable Statements .. IF (tlast) THEN ! Print solution IF (iout==2) THEN level = nlev WRITE (nout,99999) level, t WRITE (nout,99998) npts = ngpts(level) ipsol = lsol(level) DO i = 1, npts, 2 WRITE (nout,99997) xpts(i), ypts(i), sol(ipsol+i), & sol(ipsol+npts+i) END DO WRITE (nout,*) END IF END IF RETURN 99999 FORMAT (1X,'Solution at every 2nd grid point in level',I10, & ' at time ',F8.4,':') 99998 FORMAT (1X/7X,'x',10X,'y',9X,'approx c1',3X,'approx c2'/) 99997 FORMAT (2(1X,E11.4),2X,E11.4,2X,E11.4) END SUBROUTINE monit2 SUBROUTINE compute_wkspace_lens(maxlev,npde,maxpts,lenrwk,leniwk, & lenlwk) ! Returns suitable workspace lengths for the two problems ! being solved, based on trial-and-error. ! .. Implicit None Statement .. IMPLICIT NONE ! .. Scalar Arguments .. INTEGER, INTENT (OUT) :: leniwk, lenlwk, lenrwk INTEGER, INTENT (IN) :: maxlev, maxpts, npde ! .. Executable Statements .. lenrwk = 2*maxpts*npde*(5*maxlev+18*npde+9) + 2*maxpts leniwk = 2*maxpts*(14+5*maxlev) + 7*maxlev + 2 lenlwk = 2*maxpts + 400 RETURN END SUBROUTINE compute_wkspace_lens SUBROUTINE print_statistics(ts,iwk,maxlev) ! .. Implicit None Statement .. IMPLICIT NONE ! .. Scalar Arguments .. REAL (KIND=nag_wp), INTENT (IN) :: ts INTEGER, INTENT (IN) :: maxlev ! .. Array Arguments .. INTEGER, INTENT (IN) :: iwk(6*maxlev+2) ! .. Local Scalars .. INTEGER :: i, j ! .. Executable Statements .. WRITE (nout,'(1X,A)') 'Statistics:' WRITE (nout,99999) 'Time = ', ts WRITE (nout,99998) 'Total number of accepted timesteps =', iwk(1) WRITE (nout,99998) 'Total number of rejected timesteps =', iwk(2) WRITE (nout,'(1X,4(/,A))') & ' T o t a l n u m b e r o f ', & ' Residual Jacobian Newton Lin sys', & ' evals evals iters iters', ' At level ' DO j = 1, maxlev IF (iwk(j+2)/=0) THEN WRITE (nout,99997) j, (iwk(j+2+i*maxlev),i=0,3) END IF END DO WRITE (nout,'(1X,3(/,A))') & ' M a x i m u m n u m b e r o f', & ' Newton iters Lin sys iters ', ' At level ' DO j = 1, maxlev IF (iwk(j+2)/=0) THEN WRITE (nout,99996) j, (iwk(j+2+i*maxlev),i=4,5) END IF END DO WRITE (nout,*) RETURN 99999 FORMAT (1X,A,F8.4) 99998 FORMAT (1X,A,I5) 99997 FORMAT (I8,4I10) 99996 FORMAT (I8,2I14) END SUBROUTINE print_statistics END MODULE d03rafe_mod PROGRAM d03rafe ! D03RAF Example Main Program ! .. Use Statements .. USE d03rafe_mod, ONLY : nout ! .. Implicit None Statement .. IMPLICIT NONE ! .. Executable Statements .. WRITE (nout,*) 'D03RAF Example Program Results' CALL ex1 CALL ex2 CONTAINS SUBROUTINE ex1 ! .. Use Statements .. USE nag_library, ONLY : d03raf USE d03rafe_mod, ONLY : activ_energy, bndry1, compute_wkspace_lens, & damkohler, heat_release, iout, itrace, & monit1, nag_wp, nin, npde1, one, pdedef1, & pdeiv1, print_statistics, reaction_rate, & xmax, xmin, ymax, ymin, zero ! .. Implicit None Statement .. IMPLICIT NONE ! .. Parameters .. INTEGER, PARAMETER :: opti(4) = 0 ! .. Local Scalars .. REAL (KIND=nag_wp) :: tols, tolt, tout, ts INTEGER :: ifail, ind, leniwk, lenlwk, & lenrwk, maxlev, npde, npts, & nx, ny ! .. Local Arrays .. REAL (KIND=nag_wp) :: dt(3), twant(2) REAL (KIND=nag_wp), ALLOCATABLE :: optr(:,:), rwk(:) INTEGER, ALLOCATABLE :: iwk(:) LOGICAL, ALLOCATABLE :: lwk(:) ! .. Intrinsic Functions .. INTRINSIC exp, max ! .. Executable Statements .. WRITE (nout,*) WRITE (nout,*) WRITE (nout,*) 'Example 1' WRITE (nout,*) ! Skip heading in data file READ (nin,*) READ (nin,*) npts npde = npde1 damkohler = reaction_rate*exp(activ_energy)/ & (heat_release*activ_energy) dt(1:3) = (/ 0.1E-2_nag_wp, zero, zero/) twant(1:2) = (/ 0.24_nag_wp, 0.25_nag_wp/) ts = zero ind = 0 nx = 21 ny = 21 tols = 0.5_nag_wp tolt = 0.01_nag_wp maxlev = max(opti(1),3) CALL compute_wkspace_lens(maxlev,npde,npts,lenrwk,leniwk,lenlwk) ALLOCATE (rwk(lenrwk),iwk(leniwk),lwk(lenlwk),optr(3,npde)) optr(1:3,1:npde) = one 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 d03raf(npde,ts,tout,dt,xmin,xmax,ymin,ymax,nx,ny,tols,tolt, & pdedef1,bndry1,pdeiv1,monit1,opti,optr,rwk,lenrwk,iwk,leniwk, & lwk,lenlwk,itrace,ind,ifail) CALL print_statistics(ts,iwk,maxlev) END DO RETURN END SUBROUTINE ex1 SUBROUTINE ex2 ! .. Use Statements .. USE nag_library, ONLY : d03raf USE d03rafe_mod, ONLY : bndry2, compute_wkspace_lens, iout, itrace, & monit2, nag_wp, nin, npde2, one, pdedef2, & pdeiv2, print_statistics, xmax, xmin, ymax, & ymin, zero ! .. Implicit None Statement .. IMPLICIT NONE ! .. Parameters .. INTEGER, PARAMETER :: opti(4) = (/ 4, 0, 0, 0/) ! .. Local Scalars .. REAL (KIND=nag_wp) :: tols, tolt, tout, ts INTEGER :: ifail, ind, leniwk, lenlwk, & lenrwk, maxlev, npde, npts, & nx, ny ! .. Local Arrays .. REAL (KIND=nag_wp) :: dt(3), twant(2) REAL (KIND=nag_wp), ALLOCATABLE :: optr(:,:), rwk(:) INTEGER, ALLOCATABLE :: iwk(:) LOGICAL, ALLOCATABLE :: lwk(:) ! .. Intrinsic Functions .. INTRINSIC max ! .. Executable Statements .. WRITE (nout,*) WRITE (nout,*) WRITE (nout,*) 'Example 2' WRITE (nout,*) READ (nin,*) READ (nin,*) npts npde = npde2 dt(1:3) = (/ 0.5E-3_nag_wp, 1.0E-6_nag_wp, zero/) twant(1:2) = (/ 0.01_nag_wp, 0.025_nag_wp/) ts = zero ind = 0 nx = 11 ny = 11 tols = 0.075_nag_wp tolt = 0.1_nag_wp maxlev = max(opti(1),3) CALL compute_wkspace_lens(maxlev,npde,npts,lenrwk,leniwk,lenlwk) ALLOCATE (rwk(lenrwk),iwk(leniwk),lwk(lenlwk),optr(3,npde)) optr(1,1:npde) = (/ 250.0_nag_wp, 1.5E6_nag_wp/) optr(2:3,1:npde) = one DO iout = 1, 2 tout = twant(iout) ifail = 0 CALL d03raf(npde,ts,tout,dt,xmin,xmax,ymin,ymax,nx,ny,tols,tolt, & pdedef2,bndry2,pdeiv2,monit2,opti,optr,rwk,lenrwk,iwk,leniwk, & lwk,lenlwk,itrace,ind,ifail) CALL print_statistics(ts,iwk,maxlev) END DO RETURN END SUBROUTINE ex2 END PROGRAM d03rafe