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

NAG AD Library Introduction
Example description
!   D03PP_T1W_F Example Program Text
!   Mark 29.0 Release. NAG Copyright 2023.

    Module d03pp_t1w_fe_mod

!     D03PP_T1W_F Example Program Module:
!            Parameters and User-defined Routines

!     .. Use Statements ..
      Use iso_c_binding, Only: c_ptr
      Use nagad_library, Only: exp, nagad_t1w_w_rtype, Operator (/),           &
                               Operator (<), Assignment (=), Operator (-),     &
                               Operator (>), Operator (>=), Operator (+),      &
                               Operator (*)
      Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Accessibility Statements ..
      Private
      Public                           :: bndary, 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(ad_handle,npde,npts,nxi,x,xi,u,nv,v,iuser,ruser)

!       .. Scalar Arguments ..
        Type (c_ptr), Intent (Inout)   :: ad_handle
        Integer, Intent (In)           :: npde, npts, nv, nxi
!       .. Array Arguments ..
        Type (nagad_t1w_w_rtype), Intent (Inout) :: ruser(*)
        Type (nagad_t1w_w_rtype), Intent (Out) :: u(npde,npts), v(nv)
        Type (nagad_t1w_w_rtype), Intent (In) :: x(npts), xi(nxi)
        Integer, Intent (Inout)        :: iuser(*)
!       .. Local Scalars ..
        Type (nagad_t1w_w_rtype)       :: a, b, c, e, t
        Integer                        :: i
!       .. 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(ad_handle,npde,t,x,u,ux,nv,v,vdot,p,q,r,ires,iuser,    &
        ruser)

!       .. Scalar Arguments ..
        Type (c_ptr), Intent (Inout)   :: ad_handle
        Type (nagad_t1w_w_rtype), Intent (In) :: t, x
        Integer, Intent (Inout)        :: ires
        Integer, Intent (In)           :: npde, nv
!       .. Array Arguments ..
        Type (nagad_t1w_w_rtype), Intent (Out) :: p(npde,npde), q(npde),       &
                                          r(npde)
        Type (nagad_t1w_w_rtype), Intent (Inout) :: ruser(*)
        Type (nagad_t1w_w_rtype), Intent (In) :: u(npde), ux(npde), v(nv),     &
                                          vdot(nv)
        Integer, Intent (Inout)        :: iuser(*)
!       .. Local Scalars ..
        Type (nagad_t1w_w_rtype)       :: 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(ad_handle,npde,t,u,ux,nv,v,vdot,ibnd,beta,gamma,ires,  &
        iuser,ruser)

!       .. Scalar Arguments ..
        Type (c_ptr), Intent (Inout)   :: ad_handle
        Type (nagad_t1w_w_rtype), Intent (In) :: t
        Integer, Intent (In)           :: ibnd, npde, nv
        Integer, Intent (Inout)        :: ires
!       .. Array Arguments ..
        Type (nagad_t1w_w_rtype), Intent (Out) :: beta(npde), gamma(npde)
        Type (nagad_t1w_w_rtype), Intent (Inout) :: ruser(*)
        Type (nagad_t1w_w_rtype), Intent (In) :: u(npde), ux(npde), v(nv),     &
                                          vdot(nv)
        Integer, Intent (Inout)        :: iuser(*)
!       .. Local Scalars ..
        Type (nagad_t1w_w_rtype)       :: a, b, c, e, ue, x
!       .. 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(ad_handle,t,npts,npde,x,u,r,fmon,iuser,ruser)

!       .. Scalar Arguments ..
        Type (c_ptr), Intent (Inout)   :: ad_handle
        Type (nagad_t1w_w_rtype), Intent (In) :: t
        Integer, Intent (In)           :: npde, npts
!       .. Array Arguments ..
        Type (nagad_t1w_w_rtype), Intent (Out) :: fmon(npts)
        Type (nagad_t1w_w_rtype), Intent (In) :: r(npde,npts), u(npde,npts),   &
                                          x(npts)
        Type (nagad_t1w_w_rtype), Intent (Inout) :: ruser(*)
        Integer, Intent (Inout)        :: iuser(*)
!       .. Local Scalars ..
        Type (nagad_t1w_w_rtype)       :: drdx, h
        Integer                        :: i, k, l
!       .. Intrinsic Procedures ..
        Intrinsic                      :: 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
          If (drdx<0.0_nag_wp) Then
            drdx = -drdx
          End If
          fmon(i) = drdx
        End Do
        fmon(npts) = fmon(npts-1)
        Return
      End Subroutine monitf
    End Module d03pp_t1w_fe_mod
    Program d03pp_t1w_fe

!     D03PP_T1W_F Example Main Program

!     .. Use Statements ..
      Use d03pp_t1w_fe_mod, Only: bndary, half, itrace, m, monitf, nin, nout,  &
                                  npde, nv, nxfix, nxi, pdedef, two, uvinit,   &
                                  zero
      Use iso_c_binding, Only: c_ptr
      Use nagad_library, Only: d03pp_t1w_f, d03pz_t1w_f, d53pc_t1w_k, min,     &
                               nagad_t1w_set_derivative, nagad_t1w_w_rtype,    &
                               x10aa_t1w_f, x10ab_t1w_f, Assignment (=),       &
                               Operator (>)
      Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Local Scalars ..
      Type (c_ptr)                     :: ad_handle
      Type (nagad_t1w_w_rtype)         :: con, dxmesh, tout, trmesh, ts,       &
                                          xratio
      Real (Kind=nag_wp)               :: dx, e, x0, xmid
      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 ..
      Type (nagad_t1w_w_rtype)         :: algopt(30), atol(1), rtol(1),        &
                                          ruser(1), rwsav(1100), xfix(nxfix),  &
                                          xi(nxi)
      Type (nagad_t1w_w_rtype), Allocatable :: u(:), uout(:,:,:), w(:), x(:),  &
                                          xout(:)
      Real (Kind=nag_wp), Allocatable  :: de(:)
      Real (Kind=nag_wp)               :: tol(2)
      Integer                          :: iuser(1), iwsav(505)
      Integer, Allocatable             :: iw(:)
      Logical                          :: lwsav(100)
      Character (80)                   :: cwsav(10)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: real
!     .. Executable Statements ..
      Write (nout,*) 'D03PP_T1W_F 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),uout(npde,intpts,itype),w(nw),x(npts),xout(intpts),    &
        iw(niw),de(npts))
      Read (nin,*) itol
      Read (nin,*) tol(1:2)
      atol(1) = tol(1)
      rtol(1) = tol(2)
      Read (nin,*) e
      ruser(1) = e

!     Create AD configuration data object
      ifail = 0
      Call x10aa_t1w_f(ad_handle,ifail)

!     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

!     Register variables to differentiate w.r.t.
      Call nagad_t1w_set_derivative(ruser(1),1.0E0_nag_wp)

!     Loop over output value of t
      ts = zero
      tout = zero
      Do it = 1, 1
        xmid = half + 0.1_nag_wp*real(it-1,kind=nag_wp)
        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 d03pp_t1w_f(ad_handle,npde,m,ts,tout,pdedef,bndary,uvinit,u,npts, &
          x,nv,d53pc_t1w_k,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,cwsav,lwsav,iwsav,rwsav,iuser,ruser,      &
          ifail)

        If (it==1) Then
          Write (nout,99997) tol(1), npts
          Write (nout,99995) nrmesh
          Write (nout,99994) 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))


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

      End Do
      Write (nout,99998) ts%value
      Write (nout,99996) iw(1), iw(2), iw(3), iw(5)

      Do i = 1, npts
        de(i) = u(i)%tangent
      End Do

      Write (nout,*)
      Write (nout,*) ' Derivatives calculated: First order tangents'
      Write (nout,*) ' Computational mode    : algorithmic'


      Write (nout,*) ' Derivatives (final time solution w.r.t. E):'
      Write (nout,*) '    x            u(t)        du(t)/dE'
      Do i = 1, npts
        Write (nout,99999) x(i)%value, u(i)%value, de(i)
      End Do
99999 Format (1X,2(1X,E12.5),1X,E10.2)

!     Remove computational data object
      ifail = 0
      Call x10ab_t1w_f(ad_handle,ifail)

99998 Format (' T = ',F6.3)
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 (2X,'Remeshing every',I3,' time steps',/)
99994 Format (2X,'E =',F8.3)
    End Program d03pp_t1w_fe