! D03PLF Example Program Text
! Mark 27.3 Release. NAG Copyright 2021.
Module d03plfe_mod
! D03PLF 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, odedef, pdedef, uvinit
! .. Parameters ..
Real (Kind=nag_wp), Parameter, Public :: el0 = 2.5_nag_wp
Real (Kind=nag_wp), Parameter, Public :: er0 = 0.25_nag_wp
Real (Kind=nag_wp), Parameter, Public :: gamma = 1.4_nag_wp
Real (Kind=nag_wp), Parameter, Public :: rl0 = 1.0_nag_wp
Real (Kind=nag_wp), Parameter, Public :: rr0 = 0.125_nag_wp
Integer, Parameter, Public :: itrace = 0, nin = 5, nout = 6, &
npde1 = 2, npde2 = 3, nv1 = 2, &
nv2 = 0, nxi1 = 2, nxi2 = 0
Logical, Parameter, Public :: print_statistics = .False.
Contains
Subroutine exact(t,u,npde,x,npts)
! Exact solution (for comparison and b.c. purposes)
! .. 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(npde,npts)
Real (Kind=nag_wp), Intent (In) :: x(*)
! .. Local Scalars ..
Real (Kind=nag_wp) :: f, g, pi, pi2, x1, x3
Integer :: i
! .. Intrinsic Procedures ..
Intrinsic :: cos, exp, sin
! .. Executable Statements ..
f = 0.0_nag_wp
pi = x01aaf(f)
pi2 = 2.0_nag_wp*pi
Do i = 1, npts
x1 = x(i) + t
x3 = x(i) - 3.0_nag_wp*t
f = exp(pi*x3)*sin(pi2*x3)
g = exp(-pi2*x1)*cos(pi2*x1)
u(1,i) = f + g
u(2,i) = f - g
End Do
Return
End Subroutine exact
Subroutine pdedef(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)
! .. Local Scalars ..
Integer :: i
! .. Executable Statements ..
c(1:npde) = 1.0_nag_wp
d(1:npde) = 0.0_nag_wp
s(1:npde) = 0.0_nag_wp
p(1:npde,1:npde) = 0.0_nag_wp
Do i = 1, npde
p(i,i) = 1.0_nag_wp
End Do
Return
End Subroutine pdedef
Subroutine bndry1(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 Scalars ..
Real (Kind=nag_wp) :: dudx
Integer :: i
! .. Local Arrays ..
Real (Kind=nag_wp) :: ue(2,1)
! .. Executable Statements ..
If (ibnd==0) Then
i = 1
Call exact(t,ue,npde,x(i),1)
g(1) = u(1,i) + u(2,i) - ue(1,1) - ue(2,1)
dudx = (u(1,i+1)-u(2,i+1)-u(1,i)+u(2,i))/(x(i+1)-x(i))
g(2) = vdot(1) - dudx
Else
i = npts
Call exact(t,ue,npde,x(i),1)
g(1) = u(1,i) - u(2,i) - ue(1,1) + ue(2,1)
dudx = (u(1,i)+u(2,i)-u(1,i-1)-u(2,i-1))/(x(i)-x(i-1))
g(2) = vdot(2) + 3.0_nag_wp*dudx
End If
Return
End Subroutine bndry1
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)
! .. Local Scalars ..
Real (Kind=nag_wp) :: tmpl, tmpr
! .. Executable Statements ..
tmpl = 3.0_nag_wp*(uleft(1)+uleft(2))
tmpr = uright(1) - uright(2)
flux(1) = 0.5_nag_wp*(tmpl-tmpr)
flux(2) = 0.5_nag_wp*(tmpl+tmpr)
Return
End Subroutine nmflx1
Subroutine odedef(npde,t,nv,v,vdot,nxi,xi,ucp,ucpx,ucpt,r,ires)
! .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In) :: t
Integer, Intent (Inout) :: ires
Integer, Intent (In) :: npde, nv, nxi
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Out) :: r(nv)
Real (Kind=nag_wp), Intent (In) :: ucp(npde,nxi), ucpt(npde,nxi), &
ucpx(npde,nxi), v(nv), vdot(nv), &
xi(nxi)
! .. Executable Statements ..
If (ires==-1) Then
r(1) = 0.0_nag_wp
r(2) = 0.0_nag_wp
Else
r(1) = v(1) - ucp(1,1) + ucp(2,1)
r(2) = v(2) - ucp(1,2) - ucp(2,2)
End If
Return
End Subroutine odedef
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)
! .. Executable Statements ..
If (ibnd==0) Then
g(1) = u(1,1) - rl0
g(2) = u(2,1)
g(3) = u(3,1) - el0
Else
g(1) = u(1,npts) - rr0
g(2) = u(2,npts)
g(3) = u(3,npts) - er0
End If
Return
End Subroutine bndry2
Subroutine nmflx2(npde,t,x,nv,v,uleft,uright,flux,ires)
! .. Use Statements ..
Use nag_library, Only: d03puf, d03pvf
! .. 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)
! .. Local Scalars ..
Integer :: ifail
Character (1) :: path, solver
! .. Executable Statements ..
ifail = 0
solver = 'R'
If (solver=='R') Then
! ROE scheme ..
Call d03puf(uleft,uright,gamma,flux,ifail)
Else
! OSHER scheme ..
path = 'P'
Call d03pvf(uleft,uright,gamma,path,flux,ifail)
End If
Return
End Subroutine nmflx2
Subroutine uvinit(npde,npts,x,u)
! .. Scalar Arguments ..
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 ..
Integer :: i
! .. Executable Statements ..
Do i = 1, npts
If (x(i)<0.5_nag_wp) Then
u(1,i) = rl0
u(2,i) = 0.0_nag_wp
u(3,i) = el0
Else If (x(i)==0.5_nag_wp) Then
u(1,i) = 0.5_nag_wp*(rl0+rr0)
u(2,i) = 0.0_nag_wp
u(3,i) = 0.5_nag_wp*(el0+er0)
Else
u(1,i) = rr0
u(2,i) = 0.0_nag_wp
u(3,i) = er0
End If
End Do
Return
End Subroutine uvinit
End Module d03plfe_mod
Program d03plfe
! D03PLF Example Main Program
! .. Use Statements ..
Use d03plfe_mod, Only: nout
Use nag_library, Only: nag_wp
! .. Implicit None Statement ..
Implicit None
! .. Executable Statements ..
Write (nout,*) 'D03PLF Example Program Results'
Call ex1
Call ex2
Contains
Subroutine ex1
! .. Use Statements ..
Use d03plfe_mod, Only: bndry1, exact, itrace, nin, nmflx1, npde1, nv1, &
nxi1, odedef, pdedef, print_statistics
Use nag_library, Only: d03plf
! .. Local Scalars ..
Real (Kind=nag_wp) :: errmax, lerr, lwgt, tout, ts
Integer :: i, ifail, ind, itask, itol, j, &
lenode, lisave, lrsave, neqn, &
nfuncs, niters, njacs, npde, npts, &
nsteps, nv, nwkres, nxi
Character (1) :: laopt, norm
! .. Local Arrays ..
Real (Kind=nag_wp) :: algopt(30), atol(1), rtol(1)
Real (Kind=nag_wp), Allocatable :: rsave(:), u(:), ue(:,:), x(:), &
xi(:)
Integer, Allocatable :: isave(:)
! .. Intrinsic Procedures ..
Intrinsic :: abs, int, max, real
! .. Executable Statements ..
Write (nout,*)
Write (nout,*)
Write (nout,*) 'Example 1'
Write (nout,*)
! Skip heading in data file
Read (nin,*)
Read (nin,*) npts
npde = npde1
nv = nv1
nxi = nxi1
neqn = npde*npts + nv
lisave = 25*neqn + 24
nwkres = npde*(2*npts+6*nxi+3*npde+26) + nxi + nv + 7*npts + 2
lenode = 11*neqn + 50
lrsave = 4*neqn + 11*neqn/2 + 1 + nwkres + lenode
lisave = lisave*4
lrsave = lrsave*4
Allocate (rsave(lrsave),u(neqn),ue(npde,npts),x(npts),xi(nxi), &
isave(lisave))
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
xi(1) = 0.0_nag_wp
xi(2) = 1.0_nag_wp
! Set initial values ..
ts = 0.0_nag_wp
Call exact(ts,u,npde,x,npts)
u(neqn-1) = u(1) - u(2)
u(neqn) = u(neqn-2) + u(neqn-3)
laopt = 'S'
ind = 0
itask = 1
algopt(1:30) = 0.0_nag_wp
! Theta integration
algopt(1) = 1.0_nag_wp
! Sparse matrix algebra parameters
algopt(29) = 0.1_nag_wp
algopt(30) = 1.1_nag_wp
tout = 0.5_nag_wp
! ifail: behaviour on error exit
! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
ifail = 0
Call d03plf(npde,ts,tout,pdedef,nmflx1,bndry1,u,npts,x,nv,odedef,nxi, &
xi,neqn,rtol,atol,itol,norm,laopt,algopt,rsave,lrsave,isave,lisave, &
itask,itrace,ind,ifail)
Write (nout,99992)
Write (nout,99991) npts
Write (nout,99990) rtol(1)
Write (nout,99989) atol(1)
! Calculate global error at t=tout : max (||u-ue||, over x)
! Get exact solution at t=tout
Call exact(tout,ue,npde,x,npts)
errmax = -1.0_nag_wp
Do i = 2, npts
lerr = 0.0_nag_wp
Do j = 1, npde
lwgt = rtol(1)*abs(ue(j,i)) + atol(1)
lerr = lerr + abs(u((i-1)*npde+j)-ue(j,i))/lwgt
End Do
lerr = lerr/real(npde,kind=nag_wp)
errmax = max(errmax,lerr)
End Do
Write (nout,99999)
Write (nout,99998) 100*int(errmax/100.0_nag_wp) + 100
! Print integration statistics (reasonably rounded)
If (print_statistics) Then
nsteps = 50*((isave(1)+25)/50)
nfuncs = 100*((isave(2)+50)/100)
njacs = 20*((isave(3)+10)/20)
niters = 100*((isave(5)+50)/100)
Write (nout,99997)
Write (nout,99996) nsteps
Write (nout,99995) nfuncs
Write (nout,99994) njacs
Write (nout,99993) niters
End If
Return
99999 Format (/,1X,'Integration Results:')
99998 Format (2X,'Global error is less than ',I3, &
' times the local error tolerance.')
99997 Format (/,1X,'Integration Statistics:')
99996 Format (2X,'Number of time steps (nearest 50) = ',I6)
99995 Format (2X,'Number of function evaluations (nearest 100) = ',I6)
99994 Format (2X,'Number of Jacobian evaluations (nearest 20) = ',I6)
99993 Format (2X,'Number of iterations (nearest 100) = ',I6)
99992 Format (/,1X,'Method Parameters:')
99991 Format (2X,'Number of mesh points used = ',I4)
99990 Format (2X,'Relative tolerance used = ',E10.3)
99989 Format (2X,'Absolute tolerance used = ',E10.3)
End Subroutine ex1
Subroutine ex2
! .. Use Statements ..
Use d03plfe_mod, Only: bndry2, el0, er0, gamma, itrace, nin, nmflx2, &
npde2, nv2, nxi2, print_statistics, rl0, rr0, &
uvinit
Use nag_library, Only: d03pek, d03plf, d03plp
! .. Local Scalars ..
Real (Kind=nag_wp) :: d, p, tout, ts, v
Integer :: i, ifail, ind, it, itask, itol, &
lenode, lisave, lrsave, mlu, neqn, &
nfuncs, niters, njacs, npde, npts, &
nskip, nsteps, nv, nwkres, nxi
Character (1) :: laopt, norm
! .. Local Arrays ..
Real (Kind=nag_wp) :: algopt(30), atol(1), rtol(1), xi(1)
Real (Kind=nag_wp), Allocatable :: rsave(:), u(:,:), x(:)
Integer, Allocatable :: isave(:)
! .. Intrinsic Procedures ..
Intrinsic :: real
! .. Executable Statements ..
Write (nout,*)
Write (nout,*)
Write (nout,*) 'Example 2'
Write (nout,*)
Read (nin,*)
Read (nin,*) npts, nskip
npde = npde2
nv = nv2
nxi = nxi2
nwkres = npde*(2*npts+3*npde+32) + 7*npts + 4
mlu = 3*npde - 1
neqn = npde*npts + nv
lenode = 9*neqn + 50
lisave = neqn + 24
lrsave = (3*mlu+1)*neqn + nwkres + lenode
Allocate (rsave(lrsave),u(npde,npts),x(npts),isave(lisave))
! Print problem parameters
Write (nout,99997)
Write (nout,99996) gamma
Write (nout,99995) el0, er0
Write (nout,99994) rl0, rr0
! Read and print method parameters
Read (nin,*) itol
Read (nin,*) norm
Read (nin,*) atol(1), rtol(1)
Write (nout,99987)
Write (nout,99986) npts
Write (nout,99985) rtol(1)
Write (nout,99984) atol(1)
! Initialize mesh
Do i = 1, npts
x(i) = real(i-1,kind=nag_wp)/real(npts-1,kind=nag_wp)
End Do
! Initial values of variables
Call uvinit(npde,npts,x,u)
xi(1) = 0.0_nag_wp
laopt = 'B'
ind = 0
itask = 1
algopt(1:30) = 0.0_nag_wp
! 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) = 0.5E-2_nag_wp
Write (nout,99999)
Write (nout,99998)
ts = 0.0_nag_wp
tout = ts
Do it = 1, 2
tout = tout + 0.1_nag_wp
! ifail: behaviour on error exit
! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
ifail = 0
Call d03plf(npde,ts,tout,d03plp,nmflx2,bndry2,u,npts,x,nv,d03pek, &
nxi,xi,neqn,rtol,atol,itol,norm,laopt,algopt,rsave,lrsave,isave, &
lisave,itask,itrace,ind,ifail)
! Calculate density, velocity and pressure ..
Do i = 1, npts, nskip
d = u(1,i)
v = u(2,i)/d
p = d*(gamma-1.0_nag_wp)*(u(3,i)/d-0.5_nag_wp*v**2)
If (i==1) Then
Write (nout,99993) ts, x(i), d, v, p
Else
Write (nout,99992) x(i), d, v, p
End If
End Do
Write (nout,*)
End Do
! Print integration statistics (reasonably rounded)
If (print_statistics) Then
nsteps = 50*((isave(1)+25)/50)
nfuncs = 50*((isave(2)+25)/50)
njacs = isave(3)
niters = isave(5)
Write (nout,99991) nsteps
Write (nout,99990) nfuncs
Write (nout,99989) njacs
Write (nout,99988) niters
End If
Return
99999 Format (/,1X,'Solution')
99998 Format (4X,'t',6X,'x',5X,'density',1X,'velocity',1X,'pressure')
99997 Format (/,' Problem Parameter and initial conditions:')
99996 Format (' gamma =',F6.3)
99995 Format (' e(x<0.5,0) =',F6.3,' e(x>0.5,0) =',F6.3)
99994 Format (' rho(x>0.5,0) =',F6.3,' rho(x>0.5,0) =',F6.3)
99993 Format (1X,F6.3,1X,F7.4,3(2X,F7.4))
99992 Format (8X,F7.3,3(2X,F7.4))
99991 Format (/,' Number of time steps (nearest 50) = ',I6)
99990 Format (' Number of function evaluations (nearest 50) = ',I6)
99989 Format (' Number of Jacobian evaluations (nearest 1) = ',I6)
99988 Format (' Number of iterations (nearest 1) = ',I6)
99987 Format (/,' Method Parameters:')
99986 Format (' Number of mesh points used = ',I4)
99985 Format (' Relative tolerance used = ',E10.3)
99984 Format (' Absolute tolerance used = ',E10.3)
End Subroutine ex2
End Program d03plfe