! D03PFF Example Program Text ! Mark 24 Release. NAG Copyright 2012. Module d03pffe_mod ! D03PFF Example Program Module: ! Parameters and User-defined Routines ! .. Use Statements .. Use nag_library, Only: nag_wp ! .. Implicit None Statement .. Implicit None ! .. Parameters .. Integer, Parameter :: nin = 5, nout = 6, npde1 = 2, & npde2 = 1 Contains Subroutine exact(t,u,npde,x,npts) ! Exact solution (for comparison and b.c. purposes) ! .. Use Statements .. Use nag_library, Only: x01aaf ! .. Parameters .. Real (Kind=nag_wp), Parameter :: half = 0.5_nag_wp Real (Kind=nag_wp), Parameter :: two = 2.0_nag_wp ! .. 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(npts) ! .. Local Scalars .. Real (Kind=nag_wp) :: pi, px1, px2, x1, x2 Integer :: i ! .. Intrinsic Procedures .. Intrinsic :: exp, sin ! .. Executable Statements .. pi = x01aaf(pi) Do i = 1, npts x1 = x(i) + t x2 = x(i) - 3.0_nag_wp*t px1 = half*sin(two*pi*x1**2) px2 = half*sin(two*pi*x2**2) u(1,i) = half*(exp(x2)+exp(x1)+px2-px1) - t*(x1+x2) u(2,i) = (exp(x2)-exp(x1)+px2+px1) + x1*x2 + 8.0_nag_wp*t**2 End Do Return End Subroutine exact Subroutine bndry1(npde,npts,t,x,u,ibnd,g,ires) ! .. Parameters .. Real (Kind=nag_wp), Parameter :: one = 1.0_nag_wp Real (Kind=nag_wp), Parameter :: two = 2.0_nag_wp ! .. Scalar Arguments .. Real (Kind=nag_wp), Intent (In) :: t Integer, Intent (In) :: ibnd, npde, npts Integer, Intent (Inout) :: ires ! .. Array Arguments .. Real (Kind=nag_wp), Intent (Out) :: g(npde) Real (Kind=nag_wp), Intent (In) :: u(npde,3), x(npts) ! .. Local Scalars .. Real (Kind=nag_wp) :: c, exu1, exu2 ! .. Local Arrays .. Real (Kind=nag_wp) :: ue(2,1) ! .. Executable Statements .. If (ibnd==0) Then Call exact(t,ue,npde,x(1),1) c = (x(2)-x(1))/(x(3)-x(2)) exu1 = (one+c)*u(1,2) - c*u(1,3) exu2 = (one+c)*u(2,2) - c*u(2,3) g(1) = two*u(1,1) + u(2,1) - two*ue(1,1) - ue(2,1) g(2) = two*u(1,1) - u(2,1) - two*exu1 + exu2 Else Call exact(t,ue,npde,x(npts),1) c = (x(npts)-x(npts-1))/(x(npts-1)-x(npts-2)) exu1 = (one+c)*u(1,2) - c*u(1,3) exu2 = (one+c)*u(2,2) - c*u(2,3) g(1) = two*u(1,1) - u(2,1) - two*ue(1,1) + ue(2,1) g(2) = two*u(1,1) + u(2,1) - two*exu1 - exu2 End If Return End Subroutine bndry1 Subroutine nmflx1(npde,t,x,uleft,uright,flux,ires) ! .. Parameters .. Real (Kind=nag_wp), Parameter :: half = 0.5_nag_wp ! .. Scalar Arguments .. Real (Kind=nag_wp), Intent (In) :: t, x Integer, Intent (Inout) :: ires Integer, Intent (In) :: npde ! .. Array Arguments .. Real (Kind=nag_wp), Intent (Out) :: flux(npde) Real (Kind=nag_wp), Intent (In) :: uleft(npde), uright(npde) ! .. Local Scalars .. Real (Kind=nag_wp) :: ltmp, rtmp ! .. Executable Statements .. ltmp = 3.0_nag_wp*uleft(1) + 1.5_nag_wp*uleft(2) rtmp = uright(1) - half*uright(2) flux(1) = half*(ltmp-rtmp) flux(2) = ltmp + rtmp Return End Subroutine nmflx1 Subroutine pdedef(npde,t,x,u,ux,p,c,d,s,ires) ! .. Scalar Arguments .. Real (Kind=nag_wp), Intent (In) :: t, x Integer, Intent (Inout) :: ires Integer, Intent (In) :: 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) ! .. Executable Statements .. p(1,1) = 1.0_nag_wp c(1) = 0.01_nag_wp d(1) = ux(1) s(1) = u(1) Return End Subroutine pdedef Subroutine bndry2(npde,npts,t,x,u,ibnd,g,ires) ! .. Scalar Arguments .. Real (Kind=nag_wp), Intent (In) :: t Integer, Intent (In) :: ibnd, npde, npts Integer, Intent (Inout) :: ires ! .. Array Arguments .. Real (Kind=nag_wp), Intent (Out) :: g(npde) Real (Kind=nag_wp), Intent (In) :: u(npde,3), x(npts) ! .. Executable Statements .. If (ibnd==0) Then g(1) = u(1,1) - 3.0_nag_wp Else g(1) = u(1,1) - 5.0_nag_wp End If Return End Subroutine bndry2 Subroutine nmflx2(npde,t,x,uleft,uright,flux,ires) ! .. Scalar Arguments .. Real (Kind=nag_wp), Intent (In) :: t, x Integer, Intent (Inout) :: ires Integer, Intent (In) :: npde ! .. Array Arguments .. Real (Kind=nag_wp), Intent (Out) :: flux(npde) Real (Kind=nag_wp), Intent (In) :: uleft(npde), uright(npde) ! .. Executable Statements .. If (x>=0.0E0_nag_wp) Then flux(1) = x*uleft(1) Else flux(1) = x*uright(1) End If Return End Subroutine nmflx2 End Module d03pffe_mod Program d03pffe ! D03PFF Example Main Program ! .. Use Statements .. Use d03pffe_mod, Only: nout ! .. Implicit None Statement .. Implicit None ! .. Executable Statements .. Write (nout,*) 'D03PFF Example Program Results' Call ex1 Call ex2 Contains Subroutine ex1 ! .. Use Statements .. Use nag_library, Only: d03pff, d03pfp, nag_wp Use d03pffe_mod, Only: bndry1, exact, nin, nmflx1, npde1 ! .. Local Scalars .. Real (Kind=nag_wp) :: dx, tout, ts, tsmax Integer :: i, ifail, inc, ind, it, itask, & itrace, lisave, lrsave, nop, & npde, npts, outpts ! .. Local Arrays .. Real (Kind=nag_wp) :: acc(2) Real (Kind=nag_wp), Allocatable :: rsave(:), u(:,:), ue(:,:), & x(:), 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 npde = npde1 Read (nin,*) Read (nin,*) npts, inc, outpts lisave = 24 + npde*npts lrsave = (11+9*npde)*npde*npts + (32+3*npde)*npde + 7*npts + 54 Allocate (rsave(lrsave),u(npde,npts),ue(npde,outpts),x(npts), & xout(outpts),isave(lisave)) Read (nin,*) acc(1:2) Read (nin,*) itrace Read (nin,*) tsmax ! Initialise mesh dx = 1.0_nag_wp/real(npts-1,kind=nag_wp) Do i = 1, npts x(i) = real(i-1,kind=nag_wp)*dx End Do ! Set initial values Read (nin,*) ts Call exact(ts,u,npde,x,npts) ind = 0 itask = 1 Do it = 1, 2 tout = 0.1E0_nag_wp*real(it,kind=nag_wp) ! ifail: behaviour on error exit ! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft ifail = 0 Call d03pff(npde,ts,tout,d03pfp,nmflx1,bndry1,u,npts,x,acc,tsmax, & rsave,lrsave,isave,lisave,itask,itrace,ind,ifail) If (it==1) Then Write (nout,99996) npts, acc(1), acc(2) Write (nout,99999) End If ! Set output points nop = 0 Do i = 1, npts, inc nop = nop + 1 xout(nop) = x(i) End Do Write (nout,99995) ts ! Check against exact solution Call exact(tout,ue,npde,xout,nop) nop = 0 Do i = 1, npts, inc nop = nop + 1 Write (nout,99998) xout(nop), u(1,i), ue(1,nop), u(2,i), ue(2,nop) End Do End Do Write (nout,99997) isave(1), isave(2), isave(3), isave(5) Return 99999 Format (8X,'X',8X,'Approx U',4X,'Exact U',5X,'Approx V',4X,'Exac', & 't V') 99998 Format (5(3X,F9.4)) 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 (/' NPTS = ',I4,' ACC(1) = ',E10.3,' ACC(2) = ',E10.3/) 99995 Format (/' T = ',F6.3/) End Subroutine ex1 Subroutine ex2 ! .. Use Statements .. Use nag_library, Only: d03pff, nag_wp Use d03pffe_mod, Only: bndry2, nin, nmflx2, npde2, pdedef ! .. Local Scalars .. Real (Kind=nag_wp) :: dx, tout, ts, tsmax Integer :: i, ifail, ind, it, itask, & itrace, lisave, lrsave, npde, & npts, outpts ! .. Local Arrays .. Real (Kind=nag_wp) :: acc(2) Real (Kind=nag_wp), Allocatable :: rsave(:), u(:,:), x(:), xout(:) Integer, Allocatable :: iout(:), isave(:) ! .. Intrinsic Procedures .. Intrinsic :: real ! .. Executable Statements .. Write (nout,*) Write (nout,*) Write (nout,*) 'Example 2' Write (nout,*) npde = npde2 Read (nin,*) Read (nin,*) npts, outpts lisave = 24 + npde*npts lrsave = (11+9*npde)*npde*npts + (32+3*npde)*npde + 7*npts + 54 Allocate (rsave(lrsave),u(npde,npts),x(npts),xout(outpts), & isave(lisave),iout(outpts)) Read (nin,*) iout(1:outpts) Read (nin,*) acc(1:2) Read (nin,*) itrace Read (nin,*) tsmax ! Initialise mesh dx = 2.0_nag_wp/real(npts-1,kind=nag_wp) Do i = 1, npts x(i) = -1.0_nag_wp + real(i-1,kind=nag_wp)*dx End Do ! Set initial values u(1,1:npts) = x(1:npts) + 4.0_nag_wp ind = 0 itask = 1 ! Set output points from inout indices Do i = 1, outpts xout(i) = x(iout(i)) End Do ! Loop over output value of t Read (nin,*) ts, tout Do it = 1, 2 ifail = 0 Call d03pff(npde,ts,tout,pdedef,nmflx2,bndry2,u,npts,x,acc,tsmax, & rsave,lrsave,isave,lisave,itask,itrace,ind,ifail) If (it==1) Then Write (nout,99998) npts, acc(1), acc(2) Write (nout,99996) xout(1:outpts) tout = 10.0_nag_wp End If Write (nout,99999) ts Write (nout,99995)(u(1,iout(i)),i=1,outpts) End Do Write (nout,99997) isave(1), isave(2), isave(3), isave(5) Return 99999 Format (' T = ',F6.3) 99998 Format (/' NPTS = ',I4,' ACC(1) = ',E10.3,' ACC(2) = ',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,'U ',7F9.4/) End Subroutine ex2 End Program d03pffe