! D03PSF Example Program Text
! Mark 28.3 Release. NAG Copyright 2022.
Module d03psfe_mod
! D03PSF 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, monit1, &
monit2, nmflx1, nmflx2, pdef1, &
pdef2, uvin1, uvin2
! .. Parameters ..
Real (Kind=nag_wp), Parameter :: half = 0.5_nag_wp
Real (Kind=nag_wp), Parameter, Public :: one = 1.0_nag_wp
Real (Kind=nag_wp), Parameter, Public :: zero = 0.0_nag_wp
Integer, Parameter, Public :: itrace = 0, nin = 5, nout = 6, &
npde = 1, nv = 0, nxfix = 0, nxi = 0
Contains
Subroutine exact(t,u,x,npts)
! Exact solution (for comparison and b.c. purposes)
! .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In) :: t
Integer, Intent (In) :: npts
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Out) :: u(1,npts)
Real (Kind=nag_wp), Intent (In) :: x(*)
! .. Local Scalars ..
Real (Kind=nag_wp) :: del, psi, rm, rn, s
Integer :: i
! .. Executable Statements ..
s = 0.1_nag_wp
del = 0.01_nag_wp
rm = -one/del
rn = one + s/del
Do i = 1, npts
psi = x(i) - t
If (psi<s) Then
u(1,i) = one
Else If (psi>(del+s)) Then
u(1,i) = zero
Else
u(1,i) = rm*psi + rn
End If
End Do
Return
End Subroutine exact
Subroutine uvin1(npde,npts,nxi,x,xi,u,nv,v)
! .. Use Statements ..
Use nag_library, Only: x01aaf
! .. Scalar Arguments ..
Integer, Intent (In) :: npde, npts, nv, nxi
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Out) :: u(npde,npts), v(nv)
Real (Kind=nag_wp), Intent (In) :: x(npts), xi(nxi)
! .. Local Scalars ..
Real (Kind=nag_wp) :: pi, tmp
Integer :: i
! .. Intrinsic Procedures ..
Intrinsic :: sin
! .. Executable Statements ..
tmp = zero
pi = x01aaf(tmp)
Do i = 1, npts
If (x(i)>0.2_nag_wp .And. x(i)<=0.4_nag_wp) Then
tmp = pi*(5.0_nag_wp*x(i)-one)
u(1,i) = sin(tmp)
Else
u(1,i) = zero
End If
End Do
! There are no coupled ODEs in this problem (nv = 0):
v(:) = 0._nag_wp
Return
End Subroutine uvin1
Subroutine pdef1(npde,t,x,u,ux,nv,v,vdot,p,c,d,s,ires)
! .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In) :: t, x
Integer, Intent (Inout) :: ires
Integer, Intent (In) :: npde, nv
! .. 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), v(nv), vdot(nv)
! .. Executable Statements ..
p(1,1) = one
c(1) = 0.002_nag_wp
d(1) = ux(1)
s(1) = zero
Return
End Subroutine pdef1
Subroutine bndry1(npde,npts,t,x,u,nv,v,vdot,ibnd,g,ires)
! Zero solution at both boundaries
! .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In) :: t
Integer, Intent (In) :: ibnd, npde, npts, nv
Integer, Intent (Inout) :: ires
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Out) :: g(npde)
Real (Kind=nag_wp), Intent (In) :: u(npde,npts), v(nv), vdot(nv), &
x(npts)
! .. Executable Statements ..
If (ibnd==0) Then
g(1) = u(1,1)
Else
g(1) = u(1,npts)
End If
Return
End Subroutine bndry1
Subroutine monit1(t,npts,npde,x,u,fmon)
! .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In) :: t
Integer, Intent (In) :: npde, npts
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Out) :: fmon(npts)
Real (Kind=nag_wp), Intent (In) :: u(npde,npts), x(npts)
! .. Local Scalars ..
Real (Kind=nag_wp) :: h1, h2, h3
Integer :: i
! .. Intrinsic Procedures ..
Intrinsic :: abs
! .. Executable Statements ..
Do i = 2, npts - 1
h1 = x(i) - x(i-1)
h2 = x(i+1) - x(i)
h3 = half*(x(i+1)-x(i-1))
! Second derivatives ..
fmon(i) = abs(((u(1,i+1)-u(1,i))/h2-(u(1,i)-u(1,i-1))/h1)/h3)
End Do
fmon(1) = fmon(2)
fmon(npts) = fmon(npts-1)
Return
End Subroutine monit1
Subroutine nmflx1(npde,t,x,nv,v,uleft,uright,flux,ires)
! .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In) :: t, x
Integer, Intent (Inout) :: ires
Integer, Intent (In) :: npde, nv
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Out) :: flux(npde)
Real (Kind=nag_wp), Intent (In) :: uleft(npde), uright(npde), v(nv)
! .. Executable Statements ..
flux(1) = uleft(1)
Return
End Subroutine nmflx1
Subroutine uvin2(npde,npts,nxi,x,xi,u,nv,v)
! .. Scalar Arguments ..
Integer, Intent (In) :: npde, npts, nv, nxi
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Out) :: u(npde,npts), v(nv)
Real (Kind=nag_wp), Intent (In) :: x(npts), xi(nxi)
! .. Local Scalars ..
Real (Kind=nag_wp) :: t
! .. Executable Statements ..
t = zero
Call exact(t,u,x,npts)
! There are no coupled ODEs in this problem (nv = 0):
v(:) = 0._nag_wp
Return
End Subroutine uvin2
Subroutine pdef2(npde,t,x,u,ux,nv,v,vdot,p,c,d,s,ires)
! .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In) :: t, x
Integer, Intent (Inout) :: ires
Integer, Intent (In) :: npde, nv
! .. 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), v(nv), vdot(nv)
! .. Executable Statements ..
p(1,1) = one
c(1) = zero
d(1) = zero
s(1) = -100.0_nag_wp*u(1)*(u(1)-one)*(u(1)-half)
Return
End Subroutine pdef2
Subroutine bndry2(npde,npts,t,x,u,nv,v,vdot,ibnd,g,ires)
! .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In) :: t
Integer, Intent (In) :: ibnd, npde, npts, nv
Integer, Intent (Inout) :: ires
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Out) :: g(npde)
Real (Kind=nag_wp), Intent (In) :: u(npde,npts), v(nv), vdot(nv), &
x(npts)
! .. Local Arrays ..
Real (Kind=nag_wp) :: ue(1,1)
! .. Executable Statements ..
! Solution known to be constant at both boundaries
If (ibnd==0) Then
Call exact(t,ue,x(1),1)
g(1) = ue(1,1) - u(1,1)
Else
Call exact(t,ue,x(npts),1)
g(1) = ue(1,1) - u(1,npts)
End If
Return
End Subroutine bndry2
Subroutine nmflx2(npde,t,x,nv,v,uleft,uright,flux,ires)
! .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In) :: t, x
Integer, Intent (Inout) :: ires
Integer, Intent (In) :: npde, nv
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Out) :: flux(npde)
Real (Kind=nag_wp), Intent (In) :: uleft(npde), uright(npde), v(nv)
! .. Executable Statements ..
flux(1) = uleft(1)
Return
End Subroutine nmflx2
Subroutine monit2(t,npts,npde,x,u,fmon)
! .. 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) :: fmon(npts)
Real (Kind=nag_wp), Intent (In) :: u(npde,npts), x(npts)
! .. Local Scalars ..
Real (Kind=nag_wp) :: h1, pi, ux, uxmax, xa, xl, xmax, xr, &
xx
Integer :: i
! .. Intrinsic Procedures ..
Intrinsic :: abs, cos
! .. Executable Statements ..
xx = zero
pi = x01aaf(xx)
! Locate shock ..
uxmax = zero
xmax = zero
Do i = 2, npts - 1
h1 = x(i) - x(i-1)
ux = abs((u(1,i)-u(1,i-1))/h1)
If (ux>uxmax) Then
uxmax = ux
xmax = x(i)
End If
End Do
! Desired width of nonzero region of monitor function
xa = 7.0_nag_wp/60.0_nag_wp
xl = xmax - xa
xr = xmax + xa
! Assign monitor function ..
Do i = 1, npts
If (x(i)>xl .And. x(i)<xr) Then
fmon(i) = one + cos(pi*(x(i)-xmax)/xa)
Else
fmon(i) = zero
End If
End Do
Return
End Subroutine monit2
End Module d03psfe_mod
Program d03psfe
! D03PSF Example Main Program
! .. Use Statements ..
Use d03psfe_mod, Only: nout
Use nag_library, Only: nag_wp
! .. Implicit None Statement ..
Implicit None
! .. Executable Statements ..
Write (nout,*) 'D03PSF Example Program Results'
Call ex1
Call ex2
Contains
Subroutine ex1
! .. Use Statements ..
Use d03psfe_mod, Only: bndry1, itrace, monit1, nin, nmflx1, npde, nv, &
nxfix, nxi, one, pdef1, uvin1, zero
Use nag_library, Only: d03pek, d03psf, d03pzf
! .. Local Scalars ..
Real (Kind=nag_wp) :: con, dxmesh, tout, trmesh, ts, &
xratio
Integer :: i, ifail, ind, intpts, ipminf, it, &
itask, itol, itype, lenode, lisave, &
lrsave, m, mlu, neqn, npts, nrmesh, &
nwkres
Logical :: remesh
Character (1) :: laopt, norm
! .. Local Arrays ..
Real (Kind=nag_wp) :: algopt(30), atol(1), rtol(1), &
xfix(nxfix), xi(nxi)
Real (Kind=nag_wp), Allocatable :: rsave(:), u(:,:), uout(:,:,:), &
x(:), xout(:)
Integer, Allocatable :: isave(:)
! .. Intrinsic Procedures ..
Intrinsic :: real
! .. Executable Statements ..
Write (nout,*)
Write (nout,*)
Write (nout,*) 'Example 1'
! Skip heading in data file
Read (nin,*)
Read (nin,*) npts, intpts, itype
nwkres = npde*(3*npts+3*npde+32) + 7*npts + 3
mlu = 3*npde - 1
neqn = npde*npts + nv
lenode = 11*neqn + 50
lisave = 25 + nxfix + neqn
lrsave = (3*mlu+1)*neqn + nwkres + lenode
Allocate (rsave(lrsave),u(npde,npts),uout(npde,intpts,itype),x(npts), &
xout(intpts),isave(lisave))
Read (nin,*) xout(1:intpts)
Read (nin,*) itol
Read (nin,*) norm
Read (nin,*) atol(1), rtol(1)
! Initialize mesh
Do i = 1, npts
x(i) = real(i-1,kind=nag_wp)/real(npts-1,kind=nag_wp)
End Do
! Set remesh parameters
remesh = .True.
nrmesh = 3
dxmesh = zero
trmesh = zero
con = 2.0_nag_wp/real(npts-1,kind=nag_wp)
xratio = 1.5_nag_wp
ipminf = 0
laopt = 'B'
ind = 0
itask = 1
algopt(1:30) = zero
! b.d.f. integration
algopt(1) = one
algopt(13) = 0.005_nag_wp
! Loop over output value of t
ts = zero
tout = zero
Do it = 1, 3
tout = real(it,kind=nag_wp)*0.1_nag_wp
! ifail: behaviour on error exit
! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
ifail = 0
Call d03psf(npde,ts,tout,pdef1,nmflx1,bndry1,uvin1,u,npts,x,nv, &
d03pek,nxi,xi,neqn,rtol,atol,itol,norm,laopt,algopt,remesh,nxfix, &
xfix,nrmesh,dxmesh,trmesh,ipminf,xratio,con,monit1,rsave,lrsave, &
isave,lisave,itask,itrace,ind,ifail)
If (it==1) Then
Write (nout,*)
Write (nout,99998) npts, atol, rtol
End If
Write (nout,99999) ts
Write (nout,99996) xout(1:intpts)
! Interpolate at output points ..
m = 0
ifail = 0
Call d03pzf(npde,m,u,npts,x,xout,intpts,itype,uout,ifail)
Write (nout,99995) uout(1,1:intpts,1)
End Do
Write (nout,99997) isave(1), isave(2), isave(3), isave(5)
Return
99999 Format (' T = ',F6.3)
99998 Format (/,' NPTS = ',I4,' ATOL = ',E10.3,' RTOL = ',E10.3,/)
99997 Format (' Number of integration steps in time = ',I6,/,' Number ', &
'of function evaluations = ',I6,/,' Number of Jacobian ', &
'evaluations =',I6,/,' Number of iterations = ',I6)
99996 Format (1X,'X ',7F9.4)
99995 Format (1X,'Approx U ',7F9.4,/)
End Subroutine ex1
Subroutine ex2
! .. Use Statements ..
Use d03psfe_mod, Only: bndry2, exact, itrace, monit2, nin, nmflx2, &
npde, nv, nxfix, nxi, one, pdef2, uvin2, zero
Use nag_library, Only: d03pek, d03psf, d03pzf
! .. Local Scalars ..
Real (Kind=nag_wp) :: con, dxmesh, tout, trmesh, ts, &
xratio
Integer :: i, ifail, ind, intpts, ipminf, it, &
itask, itol, itype, lenode, lisave, &
lrsave, m, mlu, neqn, npts, nrmesh, &
nwkres
Logical :: remesh
Character (1) :: laopt, norm
! .. Local Arrays ..
Real (Kind=nag_wp) :: algopt(30), atol(1), rtol(1), &
xfix(nxfix), xi(nxi)
Real (Kind=nag_wp), Allocatable :: rsave(:), u(:), ue(:,:), &
uout(:,:,:), x(:), xout(:)
Integer, Allocatable :: isave(:)
! .. Intrinsic Procedures ..
Intrinsic :: real
! .. Executable Statements ..
Write (nout,*)
Write (nout,*)
Write (nout,*) 'Example 2'
! Skip heading in data file
Read (nin,*)
Read (nin,*) npts, intpts, itype
nwkres = npde*(3*npts+3*npde+32) + 7*npts + 3
mlu = 3*npde - 1
neqn = npde*npts + nv
lenode = 11*neqn + 50
lisave = 25 + nxfix + neqn
lrsave = (3*mlu+1)*neqn + nwkres + lenode
Allocate (rsave(lrsave),u(neqn),ue(1,intpts),uout(1,intpts,itype), &
x(npts),xout(intpts),isave(lisave))
Read (nin,*) xout(1:intpts)
Read (nin,*) itol
Read (nin,*) norm
Read (nin,*) atol(1), rtol(1)
! Initialize mesh
Do i = 1, npts
x(i) = real(i-1,kind=nag_wp)/real(npts-1,kind=nag_wp)
End Do
! Set remesh parameters
remesh = .True.
nrmesh = 5
dxmesh = zero
con = one/real(npts-1,kind=nag_wp)
xratio = 1.5_nag_wp
ipminf = 0
laopt = 'B'
ind = 0
itask = 1
algopt(1:30) = zero
! Theta integration ..
algopt(1) = 2.0_nag_wp
algopt(6) = 2.0_nag_wp
algopt(7) = 2.0_nag_wp
! Max. time step ..
algopt(13) = 2.5E-3_nag_wp
ts = zero
tout = zero
Do it = 1, 2
tout = real(it,kind=nag_wp)*0.2_nag_wp
ifail = 0
Call d03psf(npde,ts,tout,pdef2,nmflx2,bndry2,uvin2,u,npts,x,nv, &
d03pek,nxi,xi,neqn,rtol,atol,itol,norm,laopt,algopt,remesh,nxfix, &
xfix,nrmesh,dxmesh,trmesh,ipminf,xratio,con,monit2,rsave,lrsave, &
isave,lisave,itask,itrace,ind,ifail)
If (it==1) Then
Write (nout,*)
Write (nout,99998) npts, atol, rtol
End If
Write (nout,99999) ts
Write (nout,99996)
! Interpolate at output points ..
m = 0
ifail = 0
Call d03pzf(npde,m,u,npts,x,xout,intpts,itype,uout,ifail)
! Check against exact solution ..
Call exact(tout,ue,xout,intpts)
Do i = 1, intpts
Write (nout,99995) xout(i), uout(1,i,1), ue(1,i)
End Do
End Do
Write (nout,99997) isave(1), isave(2), isave(3), isave(5)
Return
99999 Format (' T = ',F6.3)
99998 Format (/,' NPTS = ',I4,' ATOL = ',E10.3,' RTOL = ',E10.3,/)
99997 Format (/,' Number of integration steps in time = ',I6,/,' Number ', &
'of function evaluations = ',I6,/,' Number of Jacobian ', &
'evaluations =',I6,/,' Number of iterations = ',I6)
99996 Format (8X,'X',8X,'Approx U',4X,'Exact U',/)
99995 Format (3(3X,F9.4))
End Subroutine ex2
End Program d03psfe