! D03PLF Example Program Text ! Mark 24 Release. NAG Copyright 2012. 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 ! .. 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,ncode,v,vdot,p,c,d,s,ires) ! .. 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) ! .. 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) ! .. 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) ! .. 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) ! .. 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 ! .. 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) ! .. 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 ! .. Executable Statements .. Write (nout,*) 'D03PLF Example Program Results' Call ex1 Call ex2 Contains Subroutine ex1 ! .. Use Statements .. Use nag_library, Only: d03plf, nag_wp Use d03plfe_mod, Only: bndry1, exact, itrace, ncode1, nin, nmflx1, & npde1, nxi1, odedef, pdedef ! .. 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 Procedures .. 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, nag_wp Use d03plfe_mod, Only: bndry2, el0, er0, gamma, itrace, ncode2, nin, & nmflx2, npde2, nxi2, rl0, rr0, uvinit ! .. 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 Procedures .. 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 End Program d03plfe