! D03PLF Example Program Text
! Mark 23 Release. NAG Copyright 2011.
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
! .. Parameters ..
INTEGER, PARAMETER :: itrace = 0, ncode1 = 2, &
ncode2 = 0, nin = 5, nout = 6, &
npde1 = 2, npde2 = 3, nxi1 = 2, &
nxi2 = 0
! .. Local Scalars ..
REAL (KIND=nag_wp), SAVE :: el0, er0, gamma, rl0, rr0
CONTAINS
SUBROUTINE exact(t,u,npde,x,npts)
! Exact solution (for comparison and b.c. purposes)
! .. Use Statements ..
USE nag_library, ONLY : x01aaf
! .. Implicit None Statement ..
IMPLICIT NONE
! .. 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 Functions ..
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,ncode,v,vdot,p,c,d,s,ires)
! .. Implicit None Statement ..
IMPLICIT NONE
! .. Scalar Arguments ..
REAL (KIND=nag_wp), INTENT (IN) :: t, x
INTEGER, INTENT (INOUT) :: ires
INTEGER, INTENT (IN) :: ncode, 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), v(ncode), &
vdot(ncode)
! .. 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,ncode,v,vdot,ibnd,g,ires)
! .. Implicit None Statement ..
IMPLICIT NONE
! .. Scalar Arguments ..
REAL (KIND=nag_wp), INTENT (IN) :: t
INTEGER, INTENT (IN) :: ibnd, ncode, npde, npts
INTEGER, INTENT (INOUT) :: ires
! .. Array Arguments ..
REAL (KIND=nag_wp), INTENT (OUT) :: g(npde)
REAL (KIND=nag_wp), INTENT (IN) :: u(npde,npts), v(ncode), &
vdot(ncode), 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,ncode,v,uleft,uright,flux,ires)
! .. Implicit None Statement ..
IMPLICIT NONE
! .. Scalar Arguments ..
REAL (KIND=nag_wp), INTENT (IN) :: t, x
INTEGER, INTENT (INOUT) :: ires
INTEGER, INTENT (IN) :: ncode, npde
! .. Array Arguments ..
REAL (KIND=nag_wp), INTENT (OUT) :: flux(npde)
REAL (KIND=nag_wp), INTENT (IN) :: uleft(npde), uright(npde), &
v(ncode)
! .. 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,ncode,v,vdot,nxi,xi,ucp,ucpx,ucpt,r,ires)
! .. Implicit None Statement ..
IMPLICIT NONE
! .. Scalar Arguments ..
REAL (KIND=nag_wp), INTENT (IN) :: t
INTEGER, INTENT (INOUT) :: ires
INTEGER, INTENT (IN) :: ncode, npde, nxi
! .. Array Arguments ..
REAL (KIND=nag_wp), INTENT (OUT) :: r(ncode)
REAL (KIND=nag_wp), INTENT (IN) :: ucp(npde,*), ucpt(npde,*), &
ucpx(npde,*), v(ncode), &
vdot(ncode), 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,ncode,v,vdot,ibnd,g,ires)
! .. Implicit None Statement ..
IMPLICIT NONE
! .. Scalar Arguments ..
REAL (KIND=nag_wp), INTENT (IN) :: t
INTEGER, INTENT (IN) :: ibnd, ncode, npde, npts
INTEGER, INTENT (INOUT) :: ires
! .. Array Arguments ..
REAL (KIND=nag_wp), INTENT (OUT) :: g(npde)
REAL (KIND=nag_wp), INTENT (IN) :: u(npde,npts), v(ncode), &
vdot(ncode), 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,ncode,v,uleft,uright,flux,ires)
! .. Use Statements ..
USE nag_library, ONLY : d03puf, d03pvf
! .. Implicit None Statement ..
IMPLICIT NONE
! .. Scalar Arguments ..
REAL (KIND=nag_wp), INTENT (IN) :: t, x
INTEGER, INTENT (INOUT) :: ires
INTEGER, INTENT (IN) :: ncode, npde
! .. Array Arguments ..
REAL (KIND=nag_wp), INTENT (OUT) :: flux(npde)
REAL (KIND=nag_wp), INTENT (IN) :: uleft(npde), uright(npde), &
v(ncode)
! .. 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)
! .. Implicit None Statement ..
IMPLICIT NONE
! .. 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
! .. Implicit None Statement ..
IMPLICIT NONE
! .. External Subroutines ..
EXTERNAL ex1, ex2
! .. Executable Statements ..
WRITE (nout,*) 'D03PLF Example Program Results'
CALL ex1
CALL ex2
END PROGRAM d03plfe
SUBROUTINE ex1
! .. Use Statements ..
USE nag_library, ONLY : d03plf
USE d03plfe_mod, ONLY : bndry1, exact, itrace, nag_wp, ncode1, nin, &
nmflx1, nout, npde1, nxi1, odedef, pdedef
! .. Implicit None Statement ..
IMPLICIT NONE
! .. Local Scalars ..
REAL (KIND=nag_wp) :: tout, ts
INTEGER :: i, ifail, ind, itask, itol, &
lenode, lisave, lrsave, ncode, &
neqn, nop, npde, npts, nwkres, &
nxi, outpts
CHARACTER (1) :: laopt, norm
! .. Local Arrays ..
REAL (KIND=nag_wp) :: algopt(30), atol(1), rtol(1)
REAL (KIND=nag_wp), ALLOCATABLE :: rsave(:), u(:), ue(:,:), &
uout(:,:), x(:), xi(:), xout(:)
INTEGER, ALLOCATABLE :: isave(:)
! .. Intrinsic Functions ..
INTRINSIC 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
ncode = ncode1
nxi = nxi1
neqn = npde*npts + ncode
lisave = 25*neqn + 24
nwkres = npde*(2*npts+6*nxi+3*npde+26) + nxi + ncode + 7*npts + 2
lenode = 11*neqn + 50
lrsave = 4*neqn + 11*neqn/2 + 1 + nwkres + lenode
lisave = lisave*4
lrsave = lrsave*4
outpts = npts/20 + 1
ALLOCATE (rsave(lrsave),u(neqn),ue(npde,outpts),uout(npde,outpts), &
x(npts),xi(nxi),xout(outpts),isave(lisave))
READ (nin,*) itol
READ (nin,*) norm
READ (nin,*) atol(1), rtol(1)
! Initialise 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,ncode,odedef, &
nxi,xi,neqn,rtol,atol,itol,norm,laopt,algopt,rsave,lrsave,isave, &
lisave,itask,itrace,ind,ifail)
WRITE (nout,99995) npts, atol, rtol
! Set output points ..
nop = 0
DO i = 1, npts, 20
nop = nop + 1
xout(nop) = x(i)
END DO
WRITE (nout,99996) ts
WRITE (nout,99999)
nop = 0
DO i = 1, npts, 20
nop = nop + 1
uout(1,nop) = u(npde*(i-1)+1)
uout(2,nop) = u(npde*(i-1)+2)
END DO
! Check against exact solution ..
CALL exact(tout,ue,npde,xout,nop)
WRITE (nout,99998) (xout(i),uout(1,i),ue(1,i),uout(2,i),ue(2,i),i=1, &
nop)
WRITE (nout,99997)
WRITE (nout,99994) isave(1), isave(2), isave(3), isave(5)
RETURN
99999 FORMAT (8X,'X',8X,'Approx U1',3X,'Exact U1',4X,'Approx U2',3X, &
'Exact U2'/)
99998 FORMAT (5F12.4)
99997 FORMAT (E11.4,4E14.4)
99996 FORMAT (' T = ',F6.3)
99995 FORMAT (/' NPTS = ',I4,' ATOL = ',E10.3,' RTOL = ',E10.3/)
99994 FORMAT (' Number of integration steps in time = ',I6/' Number ', &
'of function evaluations = ',I6/' Number of Jacobian ', &
'evaluations =',I6/' Number of iterations = ',I6)
END SUBROUTINE ex1
SUBROUTINE ex2
! .. Use Statements ..
USE nag_library, ONLY : d03pek, d03plf, d03plp
USE d03plfe_mod, ONLY : bndry2, el0, er0, gamma, itrace, nag_wp, &
ncode2, nin, nmflx2, nout, npde2, nxi2, rl0, &
rr0, uvinit
! .. Implicit None Statement ..
IMPLICIT NONE
! .. Local Scalars ..
REAL (KIND=nag_wp) :: d, p, tout, ts, v
INTEGER :: i, ifail, ind, it, itask, itol, &
k, lenode, lisave, lrsave, mlu, &
ncode, neqn, npde, npts, nskip, &
nwkres, nxi, outpts
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(:,:), ue(:,:), x(:)
INTEGER, ALLOCATABLE :: isave(:)
! .. Intrinsic Functions ..
INTRINSIC real
! .. Executable Statements ..
WRITE (nout,*)
WRITE (nout,*)
WRITE (nout,*) 'Example 2'
WRITE (nout,*)
READ (nin,*)
READ (nin,*) npts, nskip
npde = npde2
ncode = ncode2
nxi = nxi2
READ (nin,*) el0, er0, gamma, rl0, rr0
nwkres = npde*(2*npts+3*npde+32) + 7*npts + 4
mlu = 3*npde - 1
neqn = npde*npts + ncode
lenode = 9*neqn + 50
lisave = neqn + 24
lrsave = (3*mlu+1)*neqn + nwkres + lenode
outpts = (npts-2*nskip-1)/nskip
ALLOCATE (rsave(lrsave),u(npde,npts),x(npts),isave(lisave), &
ue(npde,outpts))
READ (nin,*) itol
READ (nin,*) norm
READ (nin,*) atol(1), rtol(1)
WRITE (nout,99995) gamma, el0, er0, rl0, rr0
WRITE (nout,99997) npts, atol, rtol
! Initialise 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
ts = 0.0_nag_wp
WRITE (nout,99999)
DO it = 1, 2
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 d03plf(npde,ts,tout,d03plp,nmflx2,bndry2,u,npts,x,ncode,d03pek, &
nxi,xi,neqn,rtol,atol,itol,norm,laopt,algopt,rsave,lrsave,isave, &
lisave,itask,itrace,ind,ifail)
WRITE (nout,99998) ts
! Read exact data (to 4 d.p.) at output points ..
READ (nin,*)
DO i = 1, outpts
READ (nin,*) ue(1,i), ue(2,i), ue(3,i)
END DO
! Calculate density, velocity and pressure ..
k = 0
DO i = 2*nskip + 1, npts - nskip, 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)
k = k + 1
WRITE (nout,99994) x(i), d, ue(1,k), v, ue(2,k), p, ue(3,k)
END DO
END DO
WRITE (nout,99996) isave(1), isave(2), isave(3), isave(5)
RETURN
99999 FORMAT (4X,'X',5X,'APPROX D',1X,'EXACT D',2X,'APPROX V',1X,'EXAC', &
'T V',2X,'APPROX P',1X,'EXACT P')
99998 FORMAT (/' T = ',F6.3/)
99997 FORMAT (/' NPTS = ',I4,' ATOL = ',E10.3,' RTOL = ',E10.3/)
99996 FORMAT (/' Number of integration steps in time = ',I6/' Number ', &
'of function evaluations = ',I6/' Number of Jacobian ', &
'evaluations =',I6/' Number of iterations = ',I6)
99995 FORMAT (/' GAMMA =',F6.3,' EL0 =',F6.3,' ER0 =',F6.3,' RL0 =',F6.3, &
' RR0 =',F6.3)
99994 FORMAT (1X,F7.4,6(2X,F7.4))
END SUBROUTINE ex2