! D03RA_P0W_F Example Program Text
! Mark 28.3 Release. NAG Copyright 2022.
Module d03ra_p0w_fe_mod
! D03RA_P0W_F Example Program Module:
! Parameters and User-defined Routines
! .. Use Statements ..
Use iso_c_binding, Only: c_ptr
Use nag_library, Only: nag_wp
! .. 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 :: xmax = one
Real (Kind=nag_wp), Parameter, Public :: ymax = one
Real (Kind=nag_wp), Parameter, Public :: zero = 0.0_nag_wp
Real (Kind=nag_wp), Parameter, Public :: xmin = zero
Real (Kind=nag_wp), Parameter, Public :: ymin = zero
Integer, Parameter, Public :: itrace = 0, nin = 5, nout = 6, &
npde1 = 1, npde2 = 2
Contains
Subroutine pdeiv1(ad_handle,npts,npde,t,x,y,u,iuser,ruser)
! .. Scalar Arguments ..
Type (c_ptr), Intent (Inout) :: ad_handle
Real (Kind=nag_wp), Intent (In) :: t
Integer, Intent (In) :: npde, npts
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Inout) :: ruser(*)
Real (Kind=nag_wp), Intent (Out) :: u(npts,npde)
Real (Kind=nag_wp), 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)
! .. 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
! .. Scalar Arguments ..
Type (c_ptr), Intent (Inout) :: ad_handle
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 (Inout) :: ruser(*)
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)
Integer, Intent (Inout) :: iuser(*)
! .. Local Scalars ..
Real (Kind=nag_wp) :: damkohler
Integer :: i
! .. Intrinsic Procedures ..
Intrinsic :: exp
! .. Executable Statements ..
damkohler = reaction_rate*exp(activ_energy)/ &
(heat_release*activ_energy)
!$Omp Do
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
!$Omp End Do
Return
End Subroutine pdedef1
Subroutine bndry1(ad_handle,npts,npde,t,x,y,u,ut,ux,uy,nbpts,lbnd,res, &
iuser,ruser)
! .. Use Statements ..
Use nag_library, Only: x02ajf
! .. Scalar Arguments ..
Type (c_ptr), Intent (Inout) :: ad_handle
Real (Kind=nag_wp), Intent (In) :: t
Integer, Intent (In) :: nbpts, npde, npts
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Inout) :: res(npts,npde), ruser(*)
Real (Kind=nag_wp), 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
! .. Intrinsic Procedures ..
Intrinsic :: abs
! .. 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 ..
Real (Kind=nag_wp), Intent (In) :: dt, dtnew, t
Integer, Intent (Inout) :: ierr
Integer, Intent (In) :: nlev, npde
Logical, Intent (In) :: tlast
! .. Array Arguments ..
Real (Kind=nag_wp), 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
! Print solution
level = nlev
Write (nout,99999) level, t
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), ypts(k+i), sol(ipsol+i)
End Do
Write (nout,*)
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_p0w_fe_mod
Program d03ra_p0w_fe
! D03RA_P0W_F Example Main Program
! .. Use Statements ..
Use d03ra_p0w_fe_mod, Only: nout
Use nag_library, Only: nag_wp
! .. Implicit None Statement ..
Implicit None
! .. Executable Statements ..
Write (nout,*) 'D03RA_P0W_F Example Program Results'
Call ex1
Contains
Subroutine ex1
! .. Use Statements ..
Use d03ra_p0w_fe_mod, Only: bndry1, compute_wkspace_lens, itrace, &
monit1, nin, npde1, one, pdedef1, pdeiv1, &
print_statistics, xmax, xmin, ymax, ymin, &
zero
Use iso_c_binding, Only: c_ptr
Use nagad_library, Only: d03ra_p0w_f
! .. Local Scalars ..
Type (c_ptr) :: ad_handle
Real (Kind=nag_wp) :: tols, tolt, tout, ts, twant
Integer :: ifail, ind, leniwk, lenlwk, lenrwk, &
maxlev, npde, npts, nx, ny
! .. Local Arrays ..
Real (Kind=nag_wp) :: dt(3), ruser(1), rwsav(25)
Real (Kind=nag_wp), Allocatable :: optr(:,:), rwk(:)
Integer :: iuser(1), iwsav(20), opti(4)
Integer, Allocatable :: iwk(:)
Logical, Allocatable :: lwk(:)
Logical :: lwsav(5)
! .. Intrinsic Procedures ..
Intrinsic :: max
! .. Executable Statements ..
Write (nout,*)
Write (nout,*)
Write (nout,*) 'Example 1'
Write (nout,*)
! Skip heading in data file
Read (nin,*)
Read (nin,*) npts
npde = npde1
dt(1:3) = (/0.1E-2_nag_wp,zero,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
! ifail: behaviour on error exit
! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
ifail = 0
Call d03ra_p0w_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,iwk,maxlev)
Return
End Subroutine ex1
End Program d03ra_p0w_fe