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

NAG AD Library Introduction
Example description
!   D03RA_A1W_F Example Program Text
!   Mark 28.4 Release. NAG Copyright 2022.

    Module d03ra_a1w_fe_mod

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

!     .. Use Statements ..
      Use iso_c_binding, Only: c_ptr
      Use nagad_library, Only: abs, exp, nagad_a1w_w_rtype, Operator (<=),     &
                               Operator (*), Operator (+), Operator (-),       &
                               Operator (/), Assignment (=)
      Use nag_library, Only: nag_wp, x02ajf
!     .. Implicit None Statement ..
      Implicit None
!     .. Accessibility Statements ..
      Private
      Public                           :: bndry1, compute_wkspace_lens,        &
                                          monit1, pdedef1, pdeiv1,             &
                                          print_statistics
!     .. Parameters ..
      Real (Kind=nag_wp), Parameter, Public :: one = 1.0_nag_wp
      Real (Kind=nag_wp), Parameter, Public :: zero = 0.0_nag_wp
      Integer, Parameter, Public       :: itrace = 0, nin = 5, nout = 6,       &
                                          npde1 = 1, npde2 = 2
      Logical, Parameter, Public       :: lprint = .False.
    Contains
      Subroutine pdeiv1(ad_handle,npts,npde,t,x,y,u,iuser,ruser)

!       .. Scalar Arguments ..
        Type (c_ptr), Intent (Inout)   :: ad_handle
        Type (nagad_a1w_w_rtype), Intent (In) :: t
        Integer, Intent (In)           :: npde, npts
!       .. Array Arguments ..
        Type (nagad_a1w_w_rtype), Intent (Inout) :: ruser(*)
        Type (nagad_a1w_w_rtype), Intent (Out) :: u(npts,npde)
        Type (nagad_a1w_w_rtype), Intent (In) :: x(npts), y(npts)
        Integer, Intent (Inout)        :: iuser(*)
!       .. Executable Statements ..
        u(1:npts,1:npde) = one
        Return
      End Subroutine pdeiv1
      Subroutine pdedef1(ad_handle,npts,npde,t,x,y,u,ut,ux,uy,uxx,uxy,uyy,res, &
        iuser,ruser)

!       .. Scalar Arguments ..
        Type (c_ptr), Intent (Inout)   :: ad_handle
        Type (nagad_a1w_w_rtype), Intent (In) :: t
        Integer, Intent (In)           :: npde, npts
!       .. Array Arguments ..
        Type (nagad_a1w_w_rtype), Intent (Out) :: res(npts,npde)
        Type (nagad_a1w_w_rtype), Intent (Inout) :: ruser(*)
        Type (nagad_a1w_w_rtype), Intent (In) :: u(npts,npde), ut(npts,npde),  &
                                          ux(npts,npde), uxx(npts,npde),       &
                                          uxy(npts,npde), uy(npts,npde),       &
                                          uyy(npts,npde), x(npts), y(npts)
        Integer, Intent (Inout)        :: iuser(*)
!       .. Local Scalars ..
        Type (nagad_a1w_w_rtype)       :: activ_energy, damkohler, diffusion,  &
                                          heat_release, reaction_rate
        Integer                        :: i
!       .. Executable Statements ..
        activ_energy = ruser(1)
        diffusion = ruser(2)
        heat_release = ruser(3)
        reaction_rate = ruser(4)

        damkohler = reaction_rate*exp(activ_energy)/                           &
          (heat_release*activ_energy)
        Do i = 1, npts
          res(i,1:npde) = ut(i,1:npde) - diffusion*(uxx(i,1:npde)+uyy(i,1:npde &
            )) - damkohler*(1.0E0_nag_wp+heat_release-u(i,1:npde))*exp(        &
            -activ_energy/u(i,1:npde))
        End Do
        Return
      End Subroutine pdedef1
      Subroutine bndry1(ad_handle,npts,npde,t,x,y,u,ut,ux,uy,nbpts,lbnd,res,   &
        iuser,ruser)

!       .. Scalar Arguments ..
        Type (c_ptr), Intent (Inout)   :: ad_handle
        Type (nagad_a1w_w_rtype), Intent (In) :: t
        Integer, Intent (In)           :: nbpts, npde, npts
!       .. Array Arguments ..
        Type (nagad_a1w_w_rtype), Intent (Inout) :: res(npts,npde), ruser(*)
        Type (nagad_a1w_w_rtype), Intent (In) :: u(npts,npde), ut(npts,npde),  &
                                          ux(npts,npde), uy(npts,npde),        &
                                          x(npts), y(npts)
        Integer, Intent (Inout)        :: iuser(*)
        Integer, Intent (In)           :: lbnd(nbpts)
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: tol
        Integer                        :: i, j
!       .. Executable Statements ..
        tol = 10._nag_wp*x02ajf()

        Do i = 1, nbpts
          j = lbnd(i)

          If (abs(x(j))<=tol) Then
            res(j,1:npde) = ux(j,1:npde)
          Else If (abs(x(j)-one)<=tol) Then
            res(j,1:npde) = u(j,1:npde) - one
          Else If (abs(y(j))<=tol) Then
            res(j,1:npde) = uy(j,1:npde)
          Else If (abs(y(j)-one)<=tol) Then
            res(j,1:npde) = u(j,1:npde) - one
          End If
        End Do

        Return
      End Subroutine bndry1
      Subroutine monit1(npde,t,dt,dtnew,tlast,nlev,ngpts,xpts,ypts,lsol,sol,   &
        ierr)

!       .. Scalar Arguments ..
        Type (nagad_a1w_w_rtype), Intent (In) :: dt, dtnew, t
        Integer, Intent (Inout)        :: ierr
        Integer, Intent (In)           :: nlev, npde
        Logical, Intent (In)           :: tlast
!       .. Array Arguments ..
        Type (nagad_a1w_w_rtype), Intent (In) :: sol(*), xpts(*), ypts(*)
        Integer, Intent (In)           :: lsol(nlev), ngpts(nlev)
!       .. Local Scalars ..
        Integer                        :: i, ipsol, k, level, npts
!       .. Intrinsic Procedures ..
        Intrinsic                      :: sum
!       .. Executable Statements ..
        If (tlast) Then

          If (lprint) Then
!           Print solution

            level = nlev
            Write (nout,99999) level, t%value
            Write (nout,99998)
            npts = ngpts(level)
            ipsol = lsol(level)
            k = sum(ngpts(1:nlev-1))
            Do i = 1, npts, 4
              Write (nout,99997) xpts(k+i)%value, ypts(k+i)%value,             &
                sol(ipsol+i)%value
            End Do

            Write (nout,*)

          End If
        End If

        Return
99999   Format (1X,'Solution at every 4th grid point in level',I10,            &
          ' at time ',F8.4,':')
99998   Format (1X,/,7X,'x',10X,'y',8X,'approx u',/)
99997   Format (1X,1P,E11.4,2(1X,1P,E11.3))
      End Subroutine monit1
      Subroutine compute_wkspace_lens(maxlev,npde,maxpts,lenrwk,leniwk,lenlwk)

!       Returns suitable workspace lengths for the two problems
!       being solved, based on trial-and-error.

!       .. Scalar Arguments ..
        Integer, Intent (Out)          :: leniwk, lenlwk, lenrwk
        Integer, Intent (In)           :: maxlev, maxpts, npde
!       .. Executable Statements ..
        lenrwk = 2*maxpts*npde*(5*maxlev+18*npde+9) + 2*maxpts
        leniwk = 2*maxpts*(14+5*maxlev) + 7*maxlev + 2
        lenlwk = 2*maxpts + 400
        Return
      End Subroutine compute_wkspace_lens
      Subroutine print_statistics(ts,iwk,maxlev)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: ts
        Integer, Intent (In)           :: maxlev
!       .. Array Arguments ..
        Integer, Intent (In)           :: iwk(6*maxlev+2)
!       .. Local Scalars ..
        Integer                        :: i, j
!       .. Local Arrays ..
        Integer                        :: istats(4)
!       .. Executable Statements ..
        Write (nout,'(1X,A)') 'Statistics:'
        Write (nout,99999) 'Time = ', ts
        Write (nout,99998) 'Total number of accepted timesteps =', iwk(1)
        Write (nout,99998) 'Total number of rejected timesteps =', iwk(2)
        Write (nout,'(1X,4(/,A))')                                             &
          '             Total  number (rounded) of     ',                      &
          '          Residual  Jacobian    Newton   Lin sys',                  &
          '             evals     evals     iters     iters', ' At level '

        Do j = 1, maxlev
          If (iwk(j+2)/=0) Then
            istats(1:4) = iwk(j+2:j+2+3*maxlev:maxlev)
            Call round_statistics(istats)
            Write (nout,99997) j, istats(1:4)
          End If
        End Do

        Write (nout,'(1X,3(/,A))') '             Maximum  number of   ',       &
          '              Newton iters    Lin sys iters ', ' At level '

        Do j = 1, maxlev
          If (iwk(j+2)/=0) Then
            Write (nout,99996) j, (iwk(j+2+i*maxlev),i=4,5)
          End If
        End Do
        Write (nout,*)

        Return
99999   Format (1X,A,F8.4)
99998   Format (1X,A,I5)
99997   Format (I8,4I10)
99996   Format (I8,2I14)
      End Subroutine print_statistics
      Subroutine round_statistics(istat)

!       .. Array Arguments ..
        Integer, Intent (Inout)        :: istat(4)
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: lt
        Integer                        :: i, k
!       .. Intrinsic Procedures ..
        Intrinsic                      :: int, log, real
!       .. Executable Statements ..
        lt = log(10.0_nag_wp)
        Do i = 1, 4
          k = int(log(real(istat(i),kind=nag_wp))/lt)
          k = 10**k
          istat(i) = k*((istat(i)+k/2)/k)
        End Do
      End Subroutine round_statistics
    End Module d03ra_a1w_fe_mod
    Program d03ra_a1w_fe

!     D03RA_A1W_F Example Main Program

!     .. Use Statements ..
      Use d03ra_a1w_fe_mod, Only: nout
      Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Executable Statements ..
      Write (nout,*) 'D03RA_A1W_F Example Program Results'

      Call ex1

    Contains
      Subroutine ex1

!       .. Use Statements ..
        Use d03ra_a1w_fe_mod, Only: bndry1, compute_wkspace_lens, itrace,      &
                                    monit1, nin, npde1, one, pdedef1, pdeiv1,  &
                                    print_statistics, zero
        Use iso_c_binding, Only: c_ptr
        Use nagad_library, Only: d03ra_a1w_f, nagad_a1w_get_derivative,        &
                                 nagad_a1w_inc_derivative,                     &
                                 nagad_a1w_ir_interpret_adjoint_sparse,        &
                                 nagad_a1w_ir_register_variable,               &
                                 nagad_a1w_ir_remove, nagad_a1w_w_rtype,       &
                                 x10aa_a1w_f, x10ab_a1w_f, x10za_a1w_f,        &
                                 Assignment (=)
!       .. Parameters ..
        Real (Kind=nag_wp), Parameter  :: activ_energy = 20.0_nag_wp
        Real (Kind=nag_wp), Parameter  :: diffusion = 0.1_nag_wp
        Real (Kind=nag_wp), Parameter  :: heat_release = 1.0_nag_wp
        Real (Kind=nag_wp), Parameter  :: reaction_rate = 5.0_nag_wp
!       .. Local Scalars ..
        Type (c_ptr)                   :: ad_handle
        Type (nagad_a1w_w_rtype)       :: tols, tolt, tout, ts, twant, xmax,   &
                                          xmin, ymax, ymin
        Integer                        :: i, ifail, ind, ipsol, ix, iy, k,     &
                                          leniwk, lenlwk, lenrwk, lpts, lrwk,  &
                                          lsgn, lsun, maxlev, ngpts, nlev,     &
                                          npde, npts, nx, ny
!       .. Local Arrays ..
        Type (nagad_a1w_w_rtype)       :: dt(3), ruser(4), rwsav(25)
        Type (nagad_a1w_w_rtype), Allocatable :: optr(:,:), rwk(:)
        Real (Kind=nag_wp)             :: dx(4)
        Integer                        :: iuser(1), iwsav(20), opti(4)
        Integer, Allocatable           :: iwk(:)
        Logical, Allocatable           :: lwk(:)
        Logical                        :: lwsav(5)
!       .. Intrinsic Procedures ..
        Intrinsic                      :: max, sum
!       .. Executable Statements ..
        Write (nout,*)
        Write (nout,*)
        Write (nout,*) 'Example 1'
        Write (nout,*)

        ruser(1) = activ_energy
        ruser(2) = diffusion
        ruser(3) = heat_release
        ruser(4) = reaction_rate

        xmax = one
        ymax = one
        xmin = zero
        ymin = zero

!       Skip heading in data file
        Read (nin,*)
        Read (nin,*) npts

        npde = npde1

        dt(1) = 0.1E-2_nag_wp
        dt(2) = zero
        dt(3) = zero
        twant = 0.24_nag_wp
        ts = zero

!       Specify that we are starting the integration in time (ind = 0
!       normally).
!       Note: we have parallelized the loop in the function pdedef1 using
!       OpenMP so set alternative value of ind to indicate that this can be
!       run in parallel if we are using a multithreaded implementation.
!       Either option is OK for serial NAG Library implementations from
!       Mark 25 onwards.
        ind = 10

        nx = 41
        ny = 41
        tols = 0.5_nag_wp
        tolt = 0.01_nag_wp
        opti(1) = 6
        opti(2:4) = 0
        maxlev = max(opti(1),3)

        Call compute_wkspace_lens(maxlev,npde,npts,lenrwk,leniwk,lenlwk)

        Allocate (rwk(lenrwk),iwk(leniwk),lwk(lenlwk),optr(3,npde))
        optr(1:3,1:npde) = one

        tout = twant

!       Create AD tape
        Call x10za_a1w_f

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

!       Register variables to differentiate w.r.t.
        Call nagad_a1w_ir_register_variable(ruser)

!       ifail: behaviour on error exit
!              =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
        ifail = 0
        Call d03ra_a1w_f(ad_handle,npde,ts,tout,dt,xmin,xmax,ymin,ymax,nx,ny,  &
          tols,tolt,pdedef1,bndry1,pdeiv1,monit1,opti,optr,rwk,lenrwk,iwk,     &
          leniwk,lwk,lenlwk,itrace,ind,lwsav,iwsav,rwsav,iuser,ruser,ifail)

        Call print_statistics(ts%value,iwk,maxlev)

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

        ngpts = 3 + 6*maxlev
        lsgn = 4 + 8*maxlev
        lsun = 6 + 12*maxlev
        nlev = iwk(lsgn)

        Write (nout,99999) nlev, ts%value
        Write (nout,99998)
99999   Format (1X,/,1X,'Solution at mid point in level',I10,' at time ',F8.4, &
          ':')
99998   Format (1X,/,7X,'x',10X,'y',8X,'approx u',20X,'du/druser(1:4)',/)

        npts = iwk(ngpts+nlev-1)
        ipsol = iwk(lsun+nlev-1) + iwsav(2) - 1
        k = sum(iwk(ngpts:ngpts+nlev-2))
        lrwk = iwk(lsun) + 2*npts + 2
        lpts = lenrwk - lrwk + 1
        ix = k + lrwk + iwsav(2) - 1
        iy = ix + lpts/2

!       Choose mid-point
        i = npts/2
        Call nagad_a1w_inc_derivative(rwk(ipsol+i),1.0E0_nag_wp)
        Call nagad_a1w_ir_interpret_adjoint_sparse(ifail)
        dx = nagad_a1w_get_derivative(ruser)
        Write (nout,99997) rwk(ix+i-1)%value, rwk(iy+i-1)%value,               &
          rwk(ipsol+i)%value, dx

        Write (nout,*)
99997   Format (1X,1P,E11.4,6(1X,1P,E11.3))

!       Remove computational data object and tape
        ifail = 0
        Call x10ab_a1w_f(ad_handle,ifail)
        Call nagad_a1w_ir_remove

        Return
      End Subroutine ex1
    End Program d03ra_a1w_fe