! D03PFF Example Program Text
! Mark 26.1 Release. NAG Copyright 2016.
Module d03pffe_mod
! D03PFF 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, exact, nmflx1, &
nmflx2, pdedef
! .. Parameters ..
Integer, Parameter, Public :: nin = 5, nout = 6, npde1 = 2, &
npde2 = 1
Contains
Subroutine exact(t,u,npde,x,npts)
! Exact solution (for comparison and b.c. purposes)
! .. Use Statements ..
Use nag_library, Only: x01aaf
! .. Parameters ..
Real (Kind=nag_wp), Parameter :: half = 0.5_nag_wp
Real (Kind=nag_wp), Parameter :: two = 2.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) :: u(npde,npts)
Real (Kind=nag_wp), Intent (In) :: x(npts)
! .. Local Scalars ..
Real (Kind=nag_wp) :: pi, px1, px2, x1, x2
Integer :: i
! .. Intrinsic Procedures ..
Intrinsic :: exp, sin
! .. Executable Statements ..
pi = x01aaf(pi)
Do i = 1, npts
x1 = x(i) + t
x2 = x(i) - 3.0_nag_wp*t
px1 = half*sin(two*pi*x1**2)
px2 = half*sin(two*pi*x2**2)
u(1,i) = half*(exp(x2)+exp(x1)+px2-px1) - t*(x1+x2)
u(2,i) = (exp(x2)-exp(x1)+px2+px1) + x1*x2 + 8.0_nag_wp*t**2
End Do
Return
End Subroutine exact
Subroutine bndry1(npde,npts,t,x,u,ibnd,g,ires)
! .. Parameters ..
Real (Kind=nag_wp), Parameter :: one = 1.0_nag_wp
Real (Kind=nag_wp), Parameter :: two = 2.0_nag_wp
! .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In) :: t
Integer, Intent (In) :: ibnd, npde, npts
Integer, Intent (Inout) :: ires
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Out) :: g(npde)
Real (Kind=nag_wp), Intent (In) :: u(npde,3), x(npts)
! .. Local Scalars ..
Real (Kind=nag_wp) :: c, exu1, exu2
! .. Local Arrays ..
Real (Kind=nag_wp) :: ue(2,1)
! .. Executable Statements ..
If (ibnd==0) Then
Call exact(t,ue,npde,x(1),1)
c = (x(2)-x(1))/(x(3)-x(2))
exu1 = (one+c)*u(1,2) - c*u(1,3)
exu2 = (one+c)*u(2,2) - c*u(2,3)
g(1) = two*u(1,1) + u(2,1) - two*ue(1,1) - ue(2,1)
g(2) = two*u(1,1) - u(2,1) - two*exu1 + exu2
Else
Call exact(t,ue,npde,x(npts),1)
c = (x(npts)-x(npts-1))/(x(npts-1)-x(npts-2))
exu1 = (one+c)*u(1,2) - c*u(1,3)
exu2 = (one+c)*u(2,2) - c*u(2,3)
g(1) = two*u(1,1) - u(2,1) - two*ue(1,1) + ue(2,1)
g(2) = two*u(1,1) + u(2,1) - two*exu1 - exu2
End If
Return
End Subroutine bndry1
Subroutine nmflx1(npde,t,x,uleft,uright,flux,ires)
! .. Parameters ..
Real (Kind=nag_wp), Parameter :: half = 0.5_nag_wp
! .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In) :: t, x
Integer, Intent (Inout) :: ires
Integer, Intent (In) :: npde
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Out) :: flux(npde)
Real (Kind=nag_wp), Intent (In) :: uleft(npde), uright(npde)
! .. Local Scalars ..
Real (Kind=nag_wp) :: ltmp, rtmp
! .. Executable Statements ..
ltmp = 3.0_nag_wp*uleft(1) + 1.5_nag_wp*uleft(2)
rtmp = uright(1) - half*uright(2)
flux(1) = half*(ltmp-rtmp)
flux(2) = ltmp + rtmp
Return
End Subroutine nmflx1
Subroutine pdedef(npde,t,x,u,ux,p,c,d,s,ires)
! .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In) :: t, x
Integer, Intent (Inout) :: ires
Integer, Intent (In) :: npde
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Out) :: c(npde), d(npde), p(npde,npde), &
s(npde)
Real (Kind=nag_wp), Intent (In) :: u(npde), ux(npde)
! .. Executable Statements ..
p(1,1) = 1.0_nag_wp
c(1) = 0.01_nag_wp
d(1) = ux(1)
s(1) = u(1)
Return
End Subroutine pdedef
Subroutine bndry2(npde,npts,t,x,u,ibnd,g,ires)
! .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In) :: t
Integer, Intent (In) :: ibnd, npde, npts
Integer, Intent (Inout) :: ires
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Out) :: g(npde)
Real (Kind=nag_wp), Intent (In) :: u(npde,3), x(npts)
! .. Executable Statements ..
If (ibnd==0) Then
g(1) = u(1,1) - 3.0_nag_wp
Else
g(1) = u(1,1) - 5.0_nag_wp
End If
Return
End Subroutine bndry2
Subroutine nmflx2(npde,t,x,uleft,uright,flux,ires)
! .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In) :: t, x
Integer, Intent (Inout) :: ires
Integer, Intent (In) :: npde
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Out) :: flux(npde)
Real (Kind=nag_wp), Intent (In) :: uleft(npde), uright(npde)
! .. Executable Statements ..
If (x>=0.0E0_nag_wp) Then
flux(1) = x*uleft(1)
Else
flux(1) = x*uright(1)
End If
Return
End Subroutine nmflx2
End Module d03pffe_mod
Program d03pffe
! D03PFF Example Main Program
! .. Use Statements ..
Use d03pffe_mod, Only: nout
! .. Implicit None Statement ..
Implicit None
! .. Executable Statements ..
Write (nout,*) 'D03PFF Example Program Results'
Call ex1
Call ex2
Contains
Subroutine ex1
! .. Use Statements ..
Use d03pffe_mod, Only: bndry1, exact, nin, nmflx1, npde1
Use nag_library, Only: d03pff, d03pfp, nag_wp
! .. Local Scalars ..
Real (Kind=nag_wp) :: dx, tout, ts, tsmax
Integer :: i, ifail, inc, ind, it, itask, &
itrace, lisave, lrsave, nfuncs, &
niters, njacs, nop, npde, npts, &
nsteps, outpts
! .. Local Arrays ..
Real (Kind=nag_wp) :: acc(2)
Real (Kind=nag_wp), Allocatable :: rsave(:), u(:,:), ue(:,:), x(:), &
xout(:)
Integer, Allocatable :: isave(:)
! .. Intrinsic Procedures ..
Intrinsic :: real
! .. Executable Statements ..
Write (nout,*)
Write (nout,*)
Write (nout,*) 'Example 1'
Write (nout,*)
! Skip heading in data file
npde = npde1
Read (nin,*)
Read (nin,*) npts, inc, outpts
lisave = 24 + npde*npts
lrsave = (11+9*npde)*npde*npts + (32+3*npde)*npde + 7*npts + 54
Allocate (rsave(lrsave),u(npde,npts),ue(npde,outpts),x(npts), &
xout(outpts),isave(lisave))
Read (nin,*) acc(1:2)
Read (nin,*) itrace
Read (nin,*) tsmax
! Initialize mesh
dx = 1.0_nag_wp/real(npts-1,kind=nag_wp)
Do i = 1, npts
x(i) = real(i-1,kind=nag_wp)*dx
End Do
! Set initial values
Read (nin,*) ts
Call exact(ts,u,npde,x,npts)
ind = 0
itask = 1
Write (nout,99992) npts
Write (nout,99991) acc(1)
Write (nout,99990) acc(2)
Write (nout,99999)
Do it = 1, 2
tout = 0.1E0_nag_wp*real(it,kind=nag_wp)
! ifail: behaviour on error exit
! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
ifail = 0
Call d03pff(npde,ts,tout,d03pfp,nmflx1,bndry1,u,npts,x,acc,tsmax, &
rsave,lrsave,isave,lisave,itask,itrace,ind,ifail)
! Set output points
nop = 0
Do i = 1, npts, inc
nop = nop + 1
xout(nop) = x(i)
End Do
! Check against exact solution
Call exact(tout,ue,npde,xout,nop)
nop = 1
i = 1
Write (nout,99998) ts, xout(nop), u(1,i), ue(1,nop), u(2,i), &
ue(2,nop)
Do i = inc + 1, npts, inc
nop = nop + 1
Write (nout,99997) xout(nop), u(1,i), ue(1,nop), u(2,i), ue(2,nop)
End Do
End Do
! Print integration statistics (reasonably rounded)
nsteps = 5*((isave(1)+2)/5)
nfuncs = 50*((isave(2)+25)/50)
njacs = 10*((isave(3)+5)/10)
niters = 50*((isave(5)+25)/50)
Write (nout,99996) nsteps
Write (nout,99995) nfuncs
Write (nout,99994) njacs
Write (nout,99993) niters
Return
99999 Format (/,5X,'t',9X,'x',8X,'Approx u',4X,'Exact u',5X,'Approx v',4X, &
'Exact v')
99998 Format (/,1X,F6.3,5(3X,F9.4))
99997 Format (7X,5(3X,F9.4))
99996 Format (/,' Number of time steps (nearest 5) = ',I6)
99995 Format (' Number of function evaluations (nearest 50) = ',I6)
99994 Format (' Number of Jacobian evaluations (nearest 10) = ',I6)
99993 Format (' Number of iterations (nearest 50) = ',I6)
99992 Format (/,' Number of mesh points used = ',I4)
99991 Format (' Relative tolerance used = ',E10.3)
99990 Format (' Absolute tolerance used = ',E10.3)
End Subroutine ex1
Subroutine ex2
! .. Use Statements ..
Use d03pffe_mod, Only: bndry2, nin, nmflx2, npde2, pdedef
Use nag_library, Only: d03pff, nag_wp
! .. Local Scalars ..
Real (Kind=nag_wp) :: dx, tout, ts, tsmax
Integer :: i, ifail, ind, it, itask, itrace, &
lisave, lrsave, nfuncs, niters, &
njacs, npde, npts, nsteps, outpts
! .. Local Arrays ..
Real (Kind=nag_wp) :: acc(2)
Real (Kind=nag_wp), Allocatable :: rsave(:), u(:,:), x(:), xout(:)
Integer, Allocatable :: iout(:), isave(:)
! .. Intrinsic Procedures ..
Intrinsic :: real
! .. Executable Statements ..
Write (nout,*)
Write (nout,*)
Write (nout,*) 'Example 2'
Write (nout,*)
npde = npde2
Read (nin,*)
Read (nin,*) npts, outpts
lisave = 24 + npde*npts
lrsave = (11+9*npde)*npde*npts + (32+3*npde)*npde + 7*npts + 54
Allocate (rsave(lrsave),u(npde,npts),x(npts),xout(outpts), &
isave(lisave),iout(outpts))
Read (nin,*) iout(1:outpts)
Read (nin,*) acc(1:2)
Read (nin,*) itrace
Read (nin,*) tsmax
! Initialize mesh
dx = 2.0_nag_wp/real(npts-1,kind=nag_wp)
Do i = 1, npts
x(i) = -1.0_nag_wp + real(i-1,kind=nag_wp)*dx
End Do
! Set initial values
u(1,1:npts) = x(1:npts) + 4.0_nag_wp
ind = 0
itask = 1
! Set output points from inout indices
Do i = 1, outpts
xout(i) = x(iout(i))
End Do
! Two output value of t: tout read from file and 10.0
Read (nin,*) ts, tout
Write (nout,99995) npts
Write (nout,99994) acc(1)
Write (nout,99993) acc(2)
Write (nout,99992) xout(1:outpts)
Do it = 1, 2
ifail = 0
Call d03pff(npde,ts,tout,pdedef,nmflx2,bndry2,u,npts,x,acc,tsmax, &
rsave,lrsave,isave,lisave,itask,itrace,ind,ifail)
Write (nout,99991) ts, (u(1,iout(i)),i=1,outpts)
tout = 10.0_nag_wp
End Do
! Print integration statistics (reasonably rounded)
nsteps = 5*((isave(1)+2)/5)
nfuncs = 50*((isave(2)+25)/50)
njacs = 10*((isave(3)+5)/10)
niters = 50*((isave(5)+25)/50)
Write (nout,99999) nsteps
Write (nout,99998) nfuncs
Write (nout,99997) njacs
Write (nout,99996) niters
Return
99999 Format (/,' Number of time steps (nearest 5) = ',I6)
99998 Format (' Number of function evaluations (nearest 50) = ',I6)
99997 Format (' Number of Jacobian evaluations (nearest 10) = ',I6)
99996 Format (' Number of iterations (nearest 50) = ',I6)
99995 Format (/,' Number of mesh points used = ',I4)
99994 Format (' Relative tolerance used = ',E10.3)
99993 Format (' Absolute tolerance used = ',E10.3)
99992 Format (/,1X,'x',10X,7F9.4)
99991 Format (1X,'u(t=',F6.3,')',7F9.4)
End Subroutine ex2
End Program d03pffe