NAG Library Manual, Mark 28.3
Interfaces:  FL   CL   CPP   AD 

NAG FL Interface Introduction
Example description
!   D03PPA Example Program Text
!   Mark 28.3 Release. NAG Copyright 2022.

    Module d03ppae_mod

!     D03PPA 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, monitf, pdedef,       &
                                          uvinit
!     .. Parameters ..
      Real (Kind=nag_wp), Parameter    :: four = 4.0_nag_wp
      Real (Kind=nag_wp), Parameter    :: one = 1.0_nag_wp
      Real (Kind=nag_wp), Parameter    :: ptone = 0.1_nag_wp
      Real (Kind=nag_wp), Parameter, Public :: half = 0.5_nag_wp
      Real (Kind=nag_wp), Parameter, Public :: two = 2.0_nag_wp
      Real (Kind=nag_wp), Parameter, Public :: zero = 0.0_nag_wp
      Integer, Parameter, Public       :: itrace = 0, m = 0, nin = 5,          &
                                          nout = 6, npde = 1, nv = 0,          &
                                          nxfix = 0, nxi = 0
    Contains
      Subroutine uvinit(npde,npts,nxi,x,xi,u,nv,v,iuser,ruser)

!       .. Scalar Arguments ..
        Integer, Intent (In)           :: npde, npts, nv, nxi
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Inout) :: ruser(*)
        Real (Kind=nag_wp), Intent (Out) :: u(npde,npts), v(nv)
        Real (Kind=nag_wp), Intent (In) :: x(npts), xi(nxi)
        Integer, Intent (Inout)        :: iuser(*)
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: a, b, c, e, t
        Integer                        :: i
!       .. Intrinsic Procedures ..
        Intrinsic                      :: exp
!       .. Executable Statements ..
        e = ruser(1)
        t = zero
        Do i = 1, npts
          a = (x(i)-0.25_nag_wp-0.75_nag_wp*t)/(four*e)
          b = (0.9_nag_wp*x(i)-0.325_nag_wp-0.495_nag_wp*t)/(two*e)
          If (a>zero .And. a>b) Then
            a = exp(-a)
            c = (0.8_nag_wp*x(i)-0.4_nag_wp-0.24_nag_wp*t)/(four*e)
            c = exp(c)
            u(1,i) = (half+ptone*c+a)/(one+c+a)
          Else If (b>zero .And. b>=a) Then
            b = exp(-b)
            c = (-0.8_nag_wp*x(i)+0.4_nag_wp+0.24_nag_wp*t)/(four*e)
            c = exp(c)
            u(1,i) = (ptone+half*c+b)/(one+c+b)
          Else
            a = exp(a)
            b = exp(b)
            u(1,i) = (one+half*a+ptone*b)/(one+a+b)
          End If
        End Do
!       There are no coupled ODEs in this problem (nv = 0):
        v(:) = 0._nag_wp
        Return
      End Subroutine uvinit
      Subroutine pdedef(npde,t,x,u,ux,nv,v,vdot,p,q,r,ires,iuser,ruser)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: t, x
        Integer, Intent (Inout)        :: ires
        Integer, Intent (In)           :: npde, nv
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: p(npde,npde), q(npde), r(npde)
        Real (Kind=nag_wp), Intent (Inout) :: ruser(*)
        Real (Kind=nag_wp), Intent (In) :: u(npde), ux(npde), v(nv), vdot(nv)
        Integer, Intent (Inout)        :: iuser(*)
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: e
!       .. Executable Statements ..
        e = ruser(1)
        p(1,1) = one
        r(1) = e*ux(1)
        q(1) = u(1)*ux(1)
        Return
      End Subroutine pdedef
      Subroutine bndary(npde,t,u,ux,nv,v,vdot,ibnd,beta,gamma,ires,iuser,      &
        ruser)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: t
        Integer, Intent (In)           :: ibnd, npde, nv
        Integer, Intent (Inout)        :: ires
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: beta(npde), gamma(npde)
        Real (Kind=nag_wp), Intent (Inout) :: ruser(*)
        Real (Kind=nag_wp), Intent (In) :: u(npde), ux(npde), v(nv), vdot(nv)
        Integer, Intent (Inout)        :: iuser(*)
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: a, b, c, e, ue, x
!       .. Intrinsic Procedures ..
        Intrinsic                      :: exp
!       .. Executable Statements ..
        e = ruser(1)
        beta(1) = zero
        If (ibnd==0) Then
          x = zero
          a = (x-0.25_nag_wp-0.75_nag_wp*t)/(four*e)
          b = (0.9_nag_wp*x-0.325_nag_wp-0.495_nag_wp*t)/(two*e)
          If (a>zero .And. a>b) Then
            a = exp(-a)
            c = (0.8_nag_wp*x-0.4_nag_wp-0.24_nag_wp*t)/(four*e)
            c = exp(c)
            ue = (half+ptone*c+a)/(one+c+a)
          Else If (b>zero .And. b>=a) Then
            b = exp(-b)
            c = (-0.8_nag_wp*x+0.4_nag_wp+0.24_nag_wp*t)/(four*e)
            c = exp(c)
            ue = (ptone+half*c+b)/(one+c+b)
          Else
            a = exp(a)
            b = exp(b)
            ue = (one+half*a+ptone*b)/(one+a+b)
          End If
        Else
          x = one
          a = (x-0.25_nag_wp-0.75_nag_wp*t)/(four*e)
          b = (0.9_nag_wp*x-0.325_nag_wp-0.495_nag_wp*t)/(two*e)
          If (a>zero .And. a>b) Then
            a = exp(-a)
            c = (0.8_nag_wp*x-0.4_nag_wp-0.24_nag_wp*t)/(four*e)
            c = exp(c)
            ue = (half+ptone*c+a)/(one+c+a)
          Else If (b>zero .And. b>=a) Then
            b = exp(-b)
            c = (-0.8_nag_wp*x+0.4_nag_wp+0.24_nag_wp*t)/(four*e)
            c = exp(c)
            ue = (ptone+half*c+b)/(one+c+b)
          Else
            a = exp(a)
            b = exp(b)
            ue = (one+half*a+ptone*b)/(one+a+b)
          End If
        End If
        gamma(1) = u(1) - ue
        Return
      End Subroutine bndary
      Subroutine monitf(t,npts,npde,x,u,r,fmon,iuser,ruser)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: t
        Integer, Intent (In)           :: npde, npts
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: fmon(npts)
        Real (Kind=nag_wp), Intent (In) :: r(npde,npts), u(npde,npts), x(npts)
        Real (Kind=nag_wp), Intent (Inout) :: ruser(*)
        Integer, Intent (Inout)        :: iuser(*)
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: drdx, h
        Integer                        :: i, k, l
!       .. Intrinsic Procedures ..
        Intrinsic                      :: abs, max, min
!       .. Executable Statements ..
        Do i = 1, npts - 1
          k = max(1,i-1)
          l = min(npts,i+1)
          h = (x(l)-x(k))*half
!         Second derivative ..
          drdx = (r(1,i+1)-r(1,i))/h
          fmon(i) = abs(drdx)
        End Do
        fmon(npts) = fmon(npts-1)
        Return
      End Subroutine monitf
      Subroutine exact(t,x,npts,u,iuser,ruser)
!       Exact solution (for comparison purposes)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: t
        Integer, Intent (In)           :: npts
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Inout) :: ruser(*)
        Real (Kind=nag_wp), Intent (Out) :: u(npts)
        Real (Kind=nag_wp), Intent (In) :: x(npts)
        Integer, Intent (Inout)        :: iuser(*)
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: a, b, c, e
        Integer                        :: i
!       .. Intrinsic Procedures ..
        Intrinsic                      :: exp
!       .. Executable Statements ..
        e = ruser(1)
        Do i = 1, npts
          a = (x(i)-0.25_nag_wp-0.75_nag_wp*t)/(four*e)
          b = (0.9_nag_wp*x(i)-0.325_nag_wp-0.495_nag_wp*t)/(two*e)
          If (a>zero .And. a>b) Then
            a = exp(-a)
            c = (0.8_nag_wp*x(i)-0.4_nag_wp-0.24_nag_wp*t)/(four*e)
            c = exp(c)
            u(i) = (half+ptone*c+a)/(one+c+a)
          Else If (b>zero .And. b>=a) Then
            b = exp(-b)
            c = (-0.8_nag_wp*x(i)+0.4_nag_wp+0.24_nag_wp*t)/(four*e)
            c = exp(c)
            u(i) = (ptone+half*c+b)/(one+c+b)
          Else
            a = exp(a)
            b = exp(b)
            u(i) = (one+half*a+ptone*b)/(one+a+b)
          End If
        End Do
        Return
      End Subroutine exact
    End Module d03ppae_mod
    Program d03ppae

!     D03PPA Example Main Program

!     .. Use Statements ..
      Use d03ppae_mod, Only: bndary, exact, half, itrace, m, monitf, nin,      &
                             nout, npde, nv, nxfix, nxi, pdedef, two, uvinit,  &
                             zero
      Use nag_library, Only: d03ppa, d03pzf, d53pck, nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: con, dx, dxmesh, e, tout, trmesh,    &
                                          ts, x0, xmid, xratio
      Integer                          :: i, ifail, ind, intpts, ipminf, it,   &
                                          itask, itol, itype, lenode, neqn,    &
                                          niw, npts, nrmesh, nw, nwkres
      Logical                          :: remesh, theta
      Character (1)                    :: laopt, norm
!     .. Local Arrays ..
      Real (Kind=nag_wp)               :: algopt(30), atol(1), rtol(1),        &
                                          ruser(1), rwsav(1100), xfix(nxfix),  &
                                          xi(nxi)
      Real (Kind=nag_wp), Allocatable  :: u(:), ue(:), uout(:,:,:), w(:),      &
                                          x(:), xout(:)
      Integer                          :: iuser(1), iwsav(505)
      Integer, Allocatable             :: iw(:)
      Logical                          :: lwsav(100)
      Character (80)                   :: cwsav(10)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: min, real
!     .. Executable Statements ..
      Write (nout,*) 'D03PPA Example Program Results'
!     Skip heading in data file
      Read (nin,*)
      Read (nin,*) npts, intpts, itype
      neqn = npde*npts + nv
      niw = 25 + nxfix
      nwkres = npde*(npts+3*npde+21) + 7*npts + nxfix + 3
      lenode = 11*neqn + 50
      nw = neqn*neqn + neqn + nwkres + lenode

      Allocate (u(neqn),ue(intpts),uout(npde,intpts,itype),w(nw),x(npts),      &
        xout(intpts),iw(niw))
      Read (nin,*) itol
      Read (nin,*) atol(1), rtol(1)
      Read (nin,*) e
      ruser(1) = e

!     Initialise mesh
      Do i = 1, npts
        x(i) = real(i-1,kind=nag_wp)/real(npts-1,kind=nag_wp)
      End Do

!     Set remesh parameters
      remesh = .True.
      nrmesh = 3
      dxmesh = half
      con = two/real(npts-1,kind=nag_wp)
      xratio = 1.5_nag_wp
      ipminf = 0

      norm = 'A'
      laopt = 'F'
      ind = 0
      itask = 1

!     Set theta to .TRUE. if the Theta integrator is required
      theta = .False.
      algopt(1:30) = zero
      If (theta) Then
        algopt(1) = two
      End If

!     Loop over output value of t
      ts = zero
      tout = zero
      Do it = 1, 5
        xmid = half + half*tout
        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 d03ppa(npde,m,ts,tout,pdedef,bndary,uvinit,u,npts,x,nv,d53pck,    &
          nxi,xi,neqn,rtol,atol,itol,norm,laopt,algopt,remesh,nxfix,xfix,      &
          nrmesh,dxmesh,trmesh,ipminf,xratio,con,monitf,w,nw,iw,niw,itask,     &
          itrace,ind,iuser,ruser,cwsav,lwsav,iwsav,rwsav,ifail)

        If (it==1) Then
          Write (nout,99998) atol, npts
          Write (nout,99993) nrmesh
          Write (nout,99992) e
          Write (nout,*)
        End If

!       Set output points
        dx = 0.1_nag_wp
        If (tout>half) Then
          dx = 0.05_nag_wp
        End If

        x0 = xmid - half*real(intpts-1,kind=nag_wp)*dx
        Do i = 1, intpts
          xout(i) = x0
          x0 = x0 + dx
        End Do
        xout(intpts) = min(xout(intpts),x(npts))

        Write (nout,99999) ts
        Write (nout,99996) xout(1:intpts)

!       Interpolate at output points
        ifail = 0
        Call d03pzf(npde,m,u,npts,x,xout,intpts,itype,uout,ifail)

!       Check against exact solution

        Call exact(ts,xout,intpts,ue,iuser,ruser)

        Write (nout,99995) uout(1,1:intpts,1)
        Write (nout,99994) ue(1:intpts)

      End Do
      Write (nout,99997) iw(1), iw(2), iw(3), iw(5)

99999 Format (' T = ',F6.3)
99998 Format (/,/,'  Accuracy requirement =',E10.3,' Number of points = ',I3,  &
        /)
99997 Format (' Number of integration steps in time = ',I6,/,' Number o',      &
        'f function evaluations = ',I6,/,' Number of Jacobian eval',           &
        'uations =',I6,/,' Number of iterations = ',I6)
99996 Format (1X,'X           ',5F9.4)
99995 Format (1X,'Approx sol. ',5F9.4)
99994 Format (1X,'Exact  sol. ',5F9.4,/)
99993 Format (2X,'Remeshing every',I3,' time steps',/)
99992 Format (2X,'E =',F8.3)
    End Program d03ppae