! D03RAF Example Program Text ! Mark 24 Release. NAG Copyright 2012. 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) ! .. 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) ! .. 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 Procedures .. 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 ! .. 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 Procedures .. 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) ! .. 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 ! .. 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 Procedures .. 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 ! .. 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 Procedures .. 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 ! .. 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 Procedures .. 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) ! .. 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. ! .. 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) ! .. 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, nag_wp Use d03rafe_mod, Only: activ_energy, bndry1, compute_wkspace_lens, & damkohler, heat_release, iout, itrace, monit1, & nin, npde1, one, pdedef1, pdeiv1, & print_statistics, reaction_rate, xmax, xmin, & ymax, ymin, zero ! .. 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 Procedures .. 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, nag_wp Use d03rafe_mod, Only: bndry2, compute_wkspace_lens, iout, itrace, & monit2, nin, npde2, one, pdedef2, pdeiv2, & print_statistics, xmax, xmin, ymax, ymin, zero ! .. 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 Procedures .. 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