! D03RA_A1W_F Example Program Text
! Mark 28.3 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