! D03PFF Example Program Text ! Mark 23 Release. NAG Copyright 2011. 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 ! .. Implicit None Statement .. IMPLICIT NONE ! .. 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 Functions .. 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) ! .. Implicit None Statement .. IMPLICIT NONE ! .. 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) ! .. Implicit None Statement .. IMPLICIT NONE ! .. 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) ! .. Implicit None Statement .. IMPLICIT NONE ! .. 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) ! .. Implicit None Statement .. IMPLICIT NONE ! .. 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) ! .. Implicit None Statement .. IMPLICIT NONE ! .. 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 ! .. Implicit None Statement .. IMPLICIT NONE ! .. 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 Functions .. 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 ! .. Implicit None Statement .. IMPLICIT NONE ! .. 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 Functions .. 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