! D03RBF Example Program Text
! Mark 30.0 Release. NAG Copyright 2024.
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