! D03PXF Example Program Text
! Mark 25 Release. NAG Copyright 2014.
Module d03pxfe_mod
! D03PXF 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 :: bndary, numflx
! .. Parameters ..
Real (Kind=nag_wp), Parameter, Public :: alpha_l = 460.894_nag_wp
Real (Kind=nag_wp), Parameter, Public :: alpha_r = 46.095_nag_wp
Real (Kind=nag_wp), Parameter, Public :: beta_l = 19.5975_nag_wp
Real (Kind=nag_wp), Parameter, Public :: beta_r = 6.19633_nag_wp
Real (Kind=nag_wp), Parameter, Public :: half = 0.5_nag_wp
Integer, Parameter, Public :: itrace = 0, ncode = 0, nin = 5, &
nout = 6, npde = 3, nxi = 0
! .. Local Scalars ..
Real (Kind=nag_wp), Public, Save :: el0, er0, gamma, rl0, rr0, ul0, &
ur0
Contains
Subroutine bndary(npde,npts,t,x,u,ncode,v,vdot,ibnd,g,ires)
! .. 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) - ul0
g(3) = u(3,1) - el0
Else
g(1) = u(1,npts) - rr0
g(2) = u(2,npts) - ur0
g(3) = u(3,npts) - er0
End If
Return
End Subroutine bndary
Subroutine numflx(npde,t,x,ncode,v,uleft,uright,flux,ires)
! .. Use Statements ..
Use nag_library, Only: d03pxf
! .. 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) :: tol
Integer :: ifail, niter
! .. Executable Statements ..
tol = 0.0_nag_wp
niter = 0
ifail = 0
Call d03pxf(uleft,uright,gamma,tol,niter,flux,ifail)
Return
End Subroutine numflx
End Module d03pxfe_mod
Program d03pxfe
! D03PXF Example Main Program
! .. Use Statements ..
Use nag_library, Only: d03pek, d03plf, d03plp, nag_wp
Use d03pxfe_mod, Only: alpha_l, alpha_r, beta_l, beta_r, bndary, el0, &
er0, gamma, half, itrace, ncode, nin, nout, npde, &
numflx, nxi, rl0, rr0, ul0, ur0
! .. Implicit None Statement ..
Implicit None
! .. Local Scalars ..
Real (Kind=nag_wp) :: d, p, tout, ts, v
Integer :: i, ifail, ind, itask, itol, k, &
lenode, mlu, neqn, niw, npts, &
nw, nwkres
Character (1) :: laopt, norm
! .. Local Arrays ..
Real (Kind=nag_wp) :: algopt(30), atol(1), rtol(1), &
ue(3,9), xi(1)
Real (Kind=nag_wp), Allocatable :: u(:,:), w(:), x(:)
Integer, Allocatable :: iw(:)
! .. Intrinsic Procedures ..
Intrinsic :: real
! .. Executable Statements ..
Write (nout,*) 'D03PXF Example Program Results'
! Skip heading in data file
Read (nin,*)
Read (nin,*) npts
nwkres = npde*(2*npts+3*npde+32) + 7*npts + 4
mlu = 3*npde - 1
neqn = npde*npts + ncode
niw = neqn + 24
lenode = 9*neqn + 50
nw = (3*mlu+1)*neqn + nwkres + lenode
Allocate (u(npde,npts),w(nw),x(npts),iw(niw))
Read (nin,*) gamma, rl0, rr0, ul0, ur0
el0 = alpha_l/(gamma-1.0_nag_wp) + half*rl0*beta_l**2
er0 = alpha_r/(gamma-1.0_nag_wp) + half*rr0*beta_r**2
! 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
! Initial values
Do i = 1, npts
If (x(i)<half) Then
u(1,i) = rl0
u(2,i) = ul0
u(3,i) = el0
Else If (x(i)==half) Then
u(1,i) = half*(rl0+rr0)
u(2,i) = half*(ul0+ur0)
u(3,i) = half*(el0+er0)
Else
u(1,i) = rr0
u(2,i) = ur0
u(3,i) = er0
End If
End Do
Read (nin,*) itol
Read (nin,*) norm
Read (nin,*) atol(1), rtol(1)
Read (nin,*) laopt
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
tout = 0.035_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,numflx,bndary,u,npts,x,ncode,d03pek,nxi, &
xi,neqn,rtol,atol,itol,norm,laopt,algopt,w,nw,iw,niw,itask,itrace,ind, &
ifail)
Write (nout,99998) ts
Write (nout,99999)
! Read exact data at output points
Do i = 1, 9
Read (nin,*) ue(1:3,i)
End Do
! Calculate density, velocity and pressure
k = 0
Do i = 15, npts - 14, 14
d = u(1,i)
v = u(2,i)/d
p = d*(gamma-1.0_nag_wp)*(u(3,i)/d-half*v**2)
k = k + 1
Write (nout,99996) x(i), d, ue(1,k), v, ue(2,k), p, ue(3,k)
End Do
Write (nout,99997) iw(1), iw(2), iw(3), iw(5)
99999 Format (4X,'X',7X,'APPROX D',3X,'EXACT D',4X,'APPROX V',3X,'EXAC','T V', &
4X,'APPROX P',3X,'EXACT P')
99998 Format (/' T = ',F6.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,E9.2,6(E11.4))
End Program d03pxfe