! D03PEF Example Program Text
! Mark 28.3 Release. NAG Copyright 2022.
Module d03pefe_mod
! D03PEF 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, exact, pdedef, uinit
! .. Parameters ..
Real (Kind=nag_wp), Parameter :: half = 0.5_nag_wp
Real (Kind=nag_wp), Parameter :: one = 1.0_nag_wp
Integer, Parameter, Public :: nin = 5, nleft = 1, nout = 6, &
npde = 2
Contains
Subroutine uinit(npde,npts,x,u)
! Routine for PDE initial values
! .. 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
! .. Intrinsic Procedures ..
Intrinsic :: exp, sin
! .. Executable Statements ..
Do i = 1, npts
u(1,i) = exp(x(i))
u(2,i) = sin(x(i))
End Do
Return
End Subroutine uinit
Subroutine pdedef(npde,t,x,u,ut,ux,res,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) :: res(npde)
Real (Kind=nag_wp), Intent (In) :: u(npde), ut(npde), ux(npde)
! .. Executable Statements ..
If (ires==-1) Then
res(1) = ut(1)
res(2) = ut(2)
Else
res(1) = ut(1) + ux(1) + ux(2)
res(2) = ut(2) + 4.0_nag_wp*ux(1) + ux(2)
End If
Return
End Subroutine pdedef
Subroutine bndary(npde,t,ibnd,nobc,u,ut,res,ires)
! .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In) :: t
Integer, Intent (In) :: ibnd, nobc, npde
Integer, Intent (Inout) :: ires
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Out) :: res(nobc)
Real (Kind=nag_wp), Intent (In) :: u(npde), ut(npde)
! .. Local Scalars ..
Real (Kind=nag_wp) :: t1, t3
! .. Intrinsic Procedures ..
Intrinsic :: exp, sin
! .. Executable Statements ..
If (ires==-1) Then
res(1) = 0.0_nag_wp
Else If (ibnd==0) Then
t3 = -3.0_nag_wp*t
t1 = t
res(1) = u(1) - half*((exp(t3)+exp(t1))+half*(sin(t3)-sin(t1)))
Else
t3 = one - 3.0_nag_wp*t
t1 = one + t
res(1) = u(2) - ((exp(t3)-exp(t1))+half*(sin(t3)+sin(t1)))
End If
Return
End Subroutine bndary
Subroutine exact(t,npde,npts,x,u)
! Exact solution (for comparison purposes)
! .. 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) :: xt, xt3
Integer :: i
! .. Intrinsic Procedures ..
Intrinsic :: exp, sin
! .. Executable Statements ..
Do i = 1, npts
xt3 = x(i) - 3.0_nag_wp*t
xt = x(i) + t
u(1,i) = half*((exp(xt3)+exp(xt))+half*(sin(xt3)-sin(xt)))
u(2,i) = (exp(xt3)-exp(xt)) + half*(sin(xt3)+sin(xt))
End Do
Return
End Subroutine exact
End Module d03pefe_mod
Program d03pefe
! D03PEF Example Main Program
! .. Use Statements ..
Use d03pefe_mod, Only: bndary, exact, nin, nleft, nout, npde, pdedef, &
uinit
Use nag_library, Only: d03pef, nag_wp
! .. Implicit None Statement ..
Implicit None
! .. Local Scalars ..
Real (Kind=nag_wp) :: acc, tout, ts
Integer :: i, ifail, ind, it, itask, itrace, &
lisave, lrsave, neqn, npts, nwkres
! .. Local Arrays ..
Real (Kind=nag_wp), Allocatable :: eu(:,:), rsave(:), u(:,:), x(:)
Integer, Allocatable :: isave(:)
! .. Intrinsic Procedures ..
Intrinsic :: real
! .. Executable Statements ..
Write (nout,*) 'D03PEF Example Program Results'
! Skip heading in data file
Read (nin,*)
Read (nin,*) npts
nwkres = npde*(npts+21+3*npde) + 7*npts + 4
neqn = npde*npts
lisave = neqn + 24
lrsave = 11*neqn + (4*npde+nleft+2)*neqn + 50 + nwkres
Allocate (eu(npde,npts),rsave(lrsave),u(npde,npts),x(npts), &
isave(lisave))
Read (nin,*) acc
Read (nin,*) itrace
! Set spatial-mesh points
Do i = 1, npts
x(i) = real(i-1,kind=nag_wp)/real(npts-1,kind=nag_wp)
End Do
ind = 0
itask = 1
Call uinit(npde,npts,x,u)
! Loop over output value of t
Read (nin,*) ts, tout
Do it = 1, 5
tout = 0.2_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 d03pef(npde,ts,tout,pdedef,bndary,u,npts,x,nleft,acc,rsave, &
lrsave,isave,lisave,itask,itrace,ind,ifail)
If (it==1) Then
Write (nout,99997) acc, npts
Write (nout,99999) x(5), x(13), x(21), x(29), x(37)
End If
! Check against the exact solution
Call exact(tout,npde,npts,x,eu)
Write (nout,99998) ts
Write (nout,99995) u(1,5:37:8)
Write (nout,99994) eu(1,5:37:8)
Write (nout,99993) u(2,5:37:8)
Write (nout,99992) eu(2,5:37:8)
End Do
Write (nout,99996) isave(1), isave(2), isave(3), isave(5)
99999 Format (' X ',5F10.4,/)
99998 Format (' T = ',F5.2)
99997 Format (/,/,' Accuracy requirement =',E10.3,' Number of points = ',I3, &
/)
99996 Format (' Number of integration steps in time = ',I6,/,' Number o', &
'f function evaluations = ',I6,/,' Number of Jacobian eval', &
'uations =',I6,/,' Number of iterations = ',I6)
99995 Format (' Approx U1',5F10.4)
99994 Format (' Exact U1',5F10.4)
99993 Format (' Approx U2',5F10.4)
99992 Format (' Exact U2',5F10.4,/)
End Program d03pefe