! D03RAF Example Program Text
! Mark 26.2 Release. NAG Copyright 2017.
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
! .. Accessibility Statements ..
Private
Public :: bndry1, bndry2, &
compute_wkspace_lens, monit1, &
monit2, monit_dummy, pdedef1, &
pdedef2, pdeiv1, pdeiv2, &
print_statistics
! .. Parameters ..
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, Public :: one = 1.0_nag_wp
Real (Kind=nag_wp), Parameter, Public :: xmax = one
Real (Kind=nag_wp), Parameter, Public :: ymax = one
Real (Kind=nag_wp), Parameter, Public :: zero = 0.0_nag_wp
Real (Kind=nag_wp), Parameter, Public :: xmin = zero
Real (Kind=nag_wp), Parameter, Public :: ymin = zero
Integer, Parameter, Public :: itrace = 0, nin = 5, nout = 6, &
npde1 = 1, npde2 = 2
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)
! .. Parameters ..
Real (Kind=nag_wp), Parameter :: activ_energy = 20.0_nag_wp
Real (Kind=nag_wp), Parameter :: diffusion = 0.1_nag_wp
Real (Kind=nag_wp), Parameter :: heat_release = 1.0_nag_wp
Real (Kind=nag_wp), Parameter :: reaction_rate = 5.0_nag_wp
! .. 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) :: damkohler
Integer :: i
! .. Intrinsic Procedures ..
Intrinsic :: exp
! .. Executable Statements ..
damkohler = reaction_rate*exp(activ_energy)/ &
(heat_release*activ_energy)
!$Omp Do
Do i = 1, npts
res(i,1:npde) = ut(i,1:npde) - diffusion*(uxx(i,1:npde)+uyy(i,1:npde &
)) - damkohler*(1.0E0_nag_wp+heat_release-u(i,1:npde))*exp( &
-activ_energy/u(i,1:npde))
End Do
!$Omp End Do
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, k, level, npts
! .. Intrinsic Procedures ..
Intrinsic :: sum
! .. Executable Statements ..
If (tlast) Then
! Print solution
level = nlev
Write (nout,99999) level, t
Write (nout,99998)
npts = ngpts(level)
ipsol = lsol(level)
k = sum(ngpts(1:nlev-1))
Do i = 1, npts, 4
Write (nout,99997) xpts(k+i), ypts(k+i), sol(ipsol+i)
End Do
Write (nout,*)
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 (1X,1P,E11.4,2(1X,1P,E11.3))
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) :: b1, fourpi
Integer :: i
! .. Intrinsic Procedures ..
Intrinsic :: sin
! .. Executable Statements ..
fourpi = 4.0_nag_wp*x01aaf(fourpi)
Do i = 1, npts
b1 = 1.0E0_nag_wp + alpha*x(i)*y(i) + beta*sin(fourpi*x(i))*sin( &
fourpi*y(i))
res(i,1) = ut(i,1) - (uxx(i,1)+uyy(i,1)) - &
u(i,1)*(b1-u(i,1)-0.5E-6_nag_wp*u(i,2))
res(i,2) = -0.05E0_nag_wp*(uxx(i,2)+uyy(i,2)) - &
u(i,2)*(-b1+1.0E4_nag_wp*u(i,1)-u(i,2))
End Do
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, k, level, npts
! .. Intrinsic Procedures ..
Intrinsic :: sum
! .. Executable Statements ..
If (tlast) Then
! Print solution
level = nlev
Write (nout,99999) level, t
Write (nout,99998)
npts = ngpts(level)
ipsol = lsol(level)
k = sum(ngpts(1:nlev-1))
Do i = 1, npts, 2
Write (nout,99997) xpts(k+i), ypts(k+i), sol(ipsol+i), &
sol(ipsol+npts+i)
End Do
Write (nout,*)
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 (1P,2(1X,E11.3),2X,E11.3,2X,E11.3)
End Subroutine monit2
Subroutine monit_dummy(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)
! .. Executable Statements ..
Return
End Subroutine monit_dummy
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
! .. Local Arrays ..
Integer :: istats(4)
! .. 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))') &
' Total number (rounded) of ', &
' Residual Jacobian Newton Lin sys', &
' evals evals iters iters', ' At level '
Do j = 1, maxlev
If (iwk(j+2)/=0) Then
istats(1:4) = iwk(j+2:j+2+3*maxlev:maxlev)
Call round_statistics(istats)
Write (nout,99997) j, istats(1:4)
End If
End Do
Write (nout,'(1X,3(/,A))') ' Maximum number of ', &
' 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
Subroutine round_statistics(istat)
! .. Array Arguments ..
Integer, Intent (Inout) :: istat(4)
! .. Local Scalars ..
Real (Kind=nag_wp) :: lt
Integer :: i, k
! .. Intrinsic Procedures ..
Intrinsic :: int, log, real
! .. Executable Statements ..
lt = log(10.0_nag_wp)
Do i = 1, 4
k = int(log(real(istat(i),kind=nag_wp))/lt)
k = 10**k
istat(i) = k*((istat(i)+k/2)/k)
End Do
End Subroutine round_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 d03rafe_mod, Only: bndry1, compute_wkspace_lens, itrace, monit1, &
monit_dummy, nin, npde1, one, pdedef1, pdeiv1, &
print_statistics, xmax, xmin, ymax, ymin, zero
Use nag_library, Only: d03raf, nag_wp
! .. Local Scalars ..
Real (Kind=nag_wp) :: tols, tolt, tout, ts
Integer :: i, 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(:)
Integer :: opti(4)
Logical, Allocatable :: lwk(:)
! .. Intrinsic Procedures ..
Intrinsic :: 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
dt(1:3) = (/0.1E-2_nag_wp,zero,zero/)
twant(1:2) = (/0.24_nag_wp,0.25_nag_wp/)
ts = zero
! Specify that we are starting the integration in time (ind = 0
! normally).
! Note: we have parallelized the loop in the function pdedef1 using
! OpenMP so set alternative value of ind to indicate that this can be
! run in parallel if we are using a multithreaded implementation.
! Either option is OK for serial NAG Library implementations from
! Mark 25 onwards.
ind = 10
nx = 41
ny = 41
tols = 0.5_nag_wp
tolt = 0.01_nag_wp
opti(1) = 6
opti(2:4) = 0
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 i = 1, 2
tout = twant(i)
! ifail: behaviour on error exit
! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
ifail = 0
If (i==1) Then
! Use monit_dummy to avoid output first time around
Call d03raf(npde,ts,tout,dt,xmin,xmax,ymin,ymax,nx,ny,tols,tolt, &
pdedef1,bndry1,pdeiv1,monit_dummy,opti,optr,rwk,lenrwk,iwk, &
leniwk,lwk,lenlwk,itrace,ind,ifail)
Else
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)
End If
Call print_statistics(ts,iwk,maxlev)
End Do
Return
End Subroutine ex1
Subroutine ex2
! .. Use Statements ..
Use d03rafe_mod, Only: bndry2, compute_wkspace_lens, itrace, monit2, &
monit_dummy, nin, npde2, one, pdedef2, pdeiv2, &
print_statistics, xmax, xmin, ymax, ymin, zero
Use nag_library, Only: d03raf, nag_wp
! .. Parameters ..
Integer, Parameter :: opti(4) = (/4,0,0,0/)
! .. Local Scalars ..
Real (Kind=nag_wp) :: tols, tolt, tout, ts
Integer :: i, 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
! Specify that we are starting the integration in time (ind = 0
! normally).
! Note: In this second example we have not added OpenMP directives to
! parallelize the loop in the function pdedef2. Thus the alternative
! ind=10 must not be specified here, as this will not function correctly
! if a multithreaded implementation of the NAG Library is used. Adding
! OpenMP to pdedef2, that would enable ind=10 to be used here safely, is
! left as an exercise for the reader.
ind = 0
nx = 11
ny = 11
tols = 0.1_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 i = 1, 2
tout = twant(i)
ifail = 0
If (i==1) Then
! Use monit_dummy to avoid output first time around
Call d03raf(npde,ts,tout,dt,xmin,xmax,ymin,ymax,nx,ny,tols,tolt, &
pdedef2,bndry2,pdeiv2,monit_dummy,opti,optr,rwk,lenrwk,iwk, &
leniwk,lwk,lenlwk,itrace,ind,ifail)
Else
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)
End If
Call print_statistics(ts,iwk,maxlev)
End Do
Return
End Subroutine ex2
End Program d03rafe