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

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

    Module d03rbfe_mod

!     D03RBF 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, inidom, monitr, pdedef,      &
                                          pdeiv
!     .. Parameters ..
      Real (Kind=nag_wp), Parameter, Public :: one = 1.0_nag_wp
      Real (Kind=nag_wp), Parameter, Public :: zero = 0.0_nag_wp
      Real (Kind=nag_wp), Parameter, Public :: twant(2) = (/0.25_nag_wp,one/)
      Integer, Parameter, Public       :: itrace = -1, nin = 5, nout = 6,      &
                                          npde = 2, print_stats = 0
!     .. Local Scalars ..
      Logical, Public, Save            :: do_monitr
    Contains
      Subroutine pdeiv(npts,npde,t,x,y,u)

!       .. Parameters ..
        Real (Kind=nag_wp), Parameter  :: eps = 0.001_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(npts,npde)
        Real (Kind=nag_wp), Intent (In) :: x(npts), y(npts)
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: a
        Integer                        :: i
!       .. Intrinsic Procedures ..
        Intrinsic                      :: exp
!       .. Executable Statements ..
        Do i = 1, npts
          a = (4.0_nag_wp*(y(i)-x(i))-t)/(32.0_nag_wp*eps)
          If (a<=zero) Then
            u(i,1) = 0.75_nag_wp - 0.25_nag_wp/(one+exp(a))
            u(i,2) = 0.75_nag_wp + 0.25_nag_wp/(one+exp(a))
          Else
            a = -a
            u(i,1) = 0.75_nag_wp - 0.25_nag_wp*exp(a)/(one+exp(a))
            u(i,2) = 0.75_nag_wp + 0.25_nag_wp*exp(a)/(one+exp(a))
          End If
        End Do

        Return
      End Subroutine pdeiv
      Subroutine inidom(maxpts,xmin,xmax,ymin,ymax,nx,ny,npts,nrows,nbnds,     &
        nbpts,lrow,irow,icol,llbnd,ilbnd,lbnd,ierr)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: xmax, xmin, ymax, ymin
        Integer, Intent (Inout)        :: ierr
        Integer, Intent (In)           :: maxpts
        Integer, Intent (Out)          :: nbnds, nbpts, npts, nrows, nx, ny
!       .. Array Arguments ..
        Integer, Intent (Inout)        :: icol(*), ilbnd(*), irow(*), lbnd(*), &
                                          llbnd(*), lrow(*)
!       .. Local Scalars ..
        Integer                        :: i
!       .. Local Arrays ..
        Integer                        :: icold(105), ilbndd(28), irowd(11),   &
                                          lbndd(72), llbndd(28), lrowd(11)
!       .. Executable Statements ..
        icold(1:105) = (/0,1,2,0,1,2,3,4,5,6,7,8,9,10,0,1,2,3,4,5,6,7,8,9,10,  &
          0,1,2,3,4,5,6,7,8,9,10,0,1,2,3,4,5,8,9,10,0,1,2,3,4,5,6,7,8,9,10,0,  &
          1,2,3,4,5,6,7,8,9,10,0,1,2,3,4,5,6,7,8,9,10,0,1,2,3,4,5,6,7,8,0,1,2, &
          3,4,5,6,7,8,0,1,2,3,4,5,6,7,8/)
        ilbndd(1:28) = (/1,2,3,4,1,4,1,2,3,4,3,4,1,2,12,23,34,41,14,41,12,23,  &
          34,41,43,14,21,32/)
        irowd(1:11) = (/0,1,2,3,4,5,6,7,8,9,10/)

        lbndd(1:72) = (/2,4,15,26,37,46,57,68,79,88,98,99,100,101,102,103,104, &
          96,86,85,84,83,82,70,59,48,39,28,17,6,8,9,10,11,12,13,18,29,40,49,   &
          60,72,73,74,75,76,77,67,56,45,36,25,33,32,42,52,53,43,1,97,105,87,   &
          81,3,7,71,78,14,31,51,54,34/)

        llbndd(1:28) = (/1,2,11,18,19,24,31,37,42,48,53,55,56,58,59,60,61,62,  &
          63,64,65,66,67,68,69,70,71,72/)
        lrowd(1:11) = (/1,4,15,26,37,46,57,68,79,88,97/)

        nx = 11
        ny = 11

!       Check MAXPTS against rough estimate of NPTS

        npts = nx*ny
        If (maxpts<npts) Then
          ierr = -1
          Return
        End If

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

        nrows = 11
        npts = 105
        nbnds = 28
        nbpts = 72

        Do i = 1, nrows
          lrow(i) = lrowd(i)
          irow(i) = irowd(i)
        End Do

        Do i = 1, nbnds
          llbnd(i) = llbndd(i)
          ilbnd(i) = ilbndd(i)
        End Do

        Do i = 1, nbpts
          lbnd(i) = lbndd(i)
        End Do

        Do i = 1, npts
          icol(i) = icold(i)
        End Do

        Return
      End Subroutine inidom
      Subroutine pdedef(npts,npde,t,x,y,u,ut,ux,uy,uxx,uxy,uyy,res)

!       .. Parameters ..
        Real (Kind=nag_wp), Parameter  :: eps = 1E-3_nag_wp
!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: t
        Integer, Intent (In)           :: npde, npts
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: res(npts,npde)
        Real (Kind=nag_wp), 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)
!       .. Local Scalars ..
        Integer                        :: i
!       .. Executable Statements ..
        !$Omp Do
        Do i = 1, npts
          res(i,1) = -u(i,1)*ux(i,1) - u(i,2)*uy(i,1) +                        &
            eps*(uxx(i,1)+uyy(i,1))
          res(i,2) = -u(i,1)*ux(i,2) - u(i,2)*uy(i,2) +                        &
            eps*(uxx(i,2)+uyy(i,2))
          res(i,1) = ut(i,1) - res(i,1)
          res(i,2) = ut(i,2) - res(i,2)
        End Do
        !$Omp End Do

        Return
      End Subroutine pdedef
      Subroutine bndary(npts,npde,t,x,y,u,ut,ux,uy,nbnds,nbpts,llbnd,ilbnd,    &
        lbnd,res)

!       .. Parameters ..
        Real (Kind=nag_wp), Parameter  :: eps = 1E-3_nag_wp
!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: t
        Integer, Intent (In)           :: nbnds, nbpts, npde, npts
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Inout) :: res(npts,npde)
        Real (Kind=nag_wp), Intent (In) :: u(npts,npde), ut(npts,npde),        &
                                          ux(npts,npde), uy(npts,npde),        &
                                          x(npts), y(npts)
        Integer, Intent (In)           :: ilbnd(nbnds), lbnd(nbpts),           &
                                          llbnd(nbnds)
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: a
        Integer                        :: i, k
!       .. Intrinsic Procedures ..
        Intrinsic                      :: exp
!       .. Executable Statements ..
        Do k = llbnd(1), nbpts
          i = lbnd(k)
          a = (-4.0_nag_wp*x(i)+4.0_nag_wp*y(i)-t)/(32.0_nag_wp*eps)
          If (a<=zero) Then
            res(i,1) = 0.75_nag_wp - 0.25_nag_wp/(one+exp(a))
            res(i,2) = 0.75_nag_wp + 0.25_nag_wp/(one+exp(a))
          Else
            a = -a
            res(i,1) = 0.75_nag_wp - 0.25_nag_wp*exp(a)/(one+exp(a))
            res(i,2) = 0.75_nag_wp + 0.25_nag_wp*exp(a)/(one+exp(a))
          End If
          res(i,1:2) = u(i,1:2) - res(i,1:2)
        End Do

        Return
      End Subroutine bndary
      Subroutine monitr(npde,t,dt,dtnew,tlast,nlev,xmin,ymin,dxb,dyb,lgrid,    &
        istruc,lsol,sol,ierr)

!       .. Use Statements ..
        Use nag_library, Only: d03rzf
!       .. Parameters ..
        Integer, Parameter             :: maxpts = 6000, nout = 6
!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: dt, dtnew, dxb, dyb, t, xmin, ymin
        Integer, Intent (Inout)        :: ierr
        Integer, Intent (In)           :: nlev, npde
        Logical, Intent (In)           :: tlast
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (In) :: sol(*)
        Integer, Intent (In)           :: istruc(*), lgrid(*), lsol(nlev)
!       .. Local Scalars ..
        Integer                        :: i, ifail, ipsol, level, npts
!       .. Local Arrays ..
        Real (Kind=nag_wp)             :: uex(105,2), x(maxpts), y(maxpts)
!       .. Executable Statements ..
        ifail = -1
levels: Do level = 1, nlev
          If (.Not. tlast) Then
            Exit levels
          End If
          ipsol = lsol(level)

!         Get grid information

          Call d03rzf(level,nlev,xmin,ymin,dxb,dyb,lgrid,istruc,npts,x,y,      &
            maxpts,ifail)
          If (ifail/=0) Then
            ierr = 1
            Exit levels
          End If

!         Skip printing?
          If (.Not. do_monitr .Or. level/=1) Then
            Cycle levels
          End If

!         Get exact solution
          Call pdeiv(npts,npde,t,x,y,uex)
          Write (nout,*)
          Write (nout,99999) t
          Write (nout,*)
          Write (nout,99998) 'x', 'y', 'approx u', 'exact u', 'approx v',      &
            'exact v'
          Write (nout,*)
          ipsol = lsol(level)
          Do i = 1, npts, 2
            Write (nout,99997) x(i), y(i), sol(ipsol+i), uex(i,1),             &
              sol(ipsol+npts+i), uex(i,2)
          End Do
          Write (nout,*)

        End Do levels

        Return
99999   Format (' Solution at every 2nd grid point in level 1 at time ',F8.4,  &
          ':')
99998   Format (7X,A,9X,A,6X,A,2X,A,2X,A,2X,A)
99997   Format (6(1X,F9.2))
      End Subroutine monitr
    End Module d03rbfe_mod
    Program d03rbfe

!     D03RBF Example Main Program

!     .. Use Statements ..
      Use d03rbfe_mod, Only: bndary, do_monitr, inidom, itrace, monitr, nin,   &
                             nout, npde, one, pdedef, pdeiv, print_stats,      &
                             twant, zero
      Use nag_library, Only: d03rbf, nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: tols, tolt, tout, ts
      Integer                          :: i, ifail, ind, iout, j, leniwk,      &
                                          lenlwk, lenrwk, maxlev, mxlev, npts
!     .. Local Arrays ..
      Real (Kind=nag_wp)               :: dt(3)
      Real (Kind=nag_wp), Allocatable  :: optr(:,:), rwk(:)
      Integer, Allocatable             :: iwk(:)
      Integer                          :: opti(4)
      Logical, Allocatable             :: lwk(:)
!     .. Executable Statements ..
      Write (nout,*) 'D03RBF Example Program Results'
!     Skip heading in data file
      Read (nin,*)
      Read (nin,*) npts, mxlev

      leniwk = 10*npts*(5*mxlev+14) + 2 + 7*mxlev
      lenlwk = 20*npts
      lenrwk = 20*(npts*npde*(5*mxlev+9+18*npde)+2*npts)
      Allocate (rwk(lenrwk),iwk(leniwk),lwk(lenlwk),optr(3,npde))

!     Specify that we are starting the integration in time (ind = 0 normally).
!     Note: we have parallelized the loop in the function pdedef using OpenMP
!     so set the 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

      ts = zero
      dt(1) = 0.001_nag_wp
      dt(2) = 1.0E-7_nag_wp
      dt(3) = zero
      tols = 0.1_nag_wp
      tolt = 0.05_nag_wp
      opti(1) = mxlev
      maxlev = opti(1)
      opti(2:4) = 0
      optr(1:3,1:npde) = one

!     Call main routine
      Do iout = 1, 2
        do_monitr = (iout==2)
        tout = twant(iout)

!       ifail: behaviour on error exit
!              =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
        ifail = 0
        Call d03rbf(npde,ts,tout,dt,tols,tolt,inidom,pdedef,bndary,pdeiv,      &
          monitr,opti,optr,rwk,lenrwk,iwk,leniwk,lwk,lenlwk,itrace,ind,ifail)

        If (print_stats/=0) Then
!         Print statistics

          Write (nout,99999) 'Statistics:'
          Write (nout,99998) 'Time = ', ts
          Write (nout,99997) 'Total number of accepted timesteps =', iwk(1)
          Write (nout,99997) 'Total number of rejected timesteps =', iwk(2)
          Write (nout,'(1X,4(/,A,A))') '            Total  number  of     ',   &
            '  maximum  number  of     ',                                      &
            '       Residual Jacobian    Newton   ', '      Newton   ',        &
            '        evals     evals     iters   ', '       iters   ',         &
            ' Level '

          maxlev = opti(1)
          Do j = 1, maxlev
            If (iwk(j+2)/=0) Then
              Write (nout,'(I4,4I10)') j, (iwk(j+2+i*maxlev),i=0,2),           &
                iwk(j+2+4*maxlev)
            End If
          End Do
          Write (nout,*)
        End If
      End Do

99999 Format (1X,A)
99998 Format (1X,A,F8.4)
99997 Format (1X,A,I5)
    End Program d03rbfe