# NAG Library Routine Document

## 1Purpose

d03raf integrates a system of linear or nonlinear, time-dependent partial differential equations (PDEs) in two space dimensions on a rectangular domain. The method of lines is employed to reduce the PDEs to a system of ordinary differential equations (ODEs) which are solved using a backward differentiation formula (BDF) method. The resulting system of nonlinear equations is solved using a modified Newton method and a Bi-CGSTAB iterative linear solver with ILU preconditioning. Local uniform grid refinement is used to improve the accuracy of the solution. d03raf originates from the VLUGR2 package (see Blom and Verwer (1993) and Blom et al. (1996)).

## 2Specification

Fortran Interface
 Subroutine d03raf ( npde, ts, tout, dt, xmin, xmax, ymin, ymax, nx, ny, tols, tolt, opti, optr, rwk, iwk, lwk, ind,
 Integer, Intent (In) :: npde, nx, ny, opti(4), lenrwk, leniwk, lenlwk, itrace Integer, Intent (Inout) :: iwk(leniwk), ind, ifail Real (Kind=nag_wp), Intent (In) :: tout, xmin, xmax, ymin, ymax, tols, tolt, optr(3,npde) Real (Kind=nag_wp), Intent (Inout) :: ts, dt(3), rwk(lenrwk) Logical, Intent (Out) :: lwk(lenlwk) External :: pdedef, bndary, pdeiv, monitr
#include nagmk26.h
 void d03raf_ (const Integer *npde, double *ts, const double *tout, double dt[], const double *xmin, const double *xmax, const double *ymin, const double *ymax, const Integer *nx, const Integer *ny, const double *tols, const double *tolt, void (NAG_CALL *pdedef)(const Integer *npts, const Integer *npde, const double *t, const double x[], const double y[], const double u[], const double ut[], const double ux[], const double uy[], const double uxx[], const double uxy[], const double uyy[], double res[]),void (NAG_CALL *bndary)(const Integer *npts, const Integer *npde, const double *t, const double x[], const double y[], const double u[], const double ut[], const double ux[], const double uy[], const Integer *nbpts, const Integer lbnd[], double res[]),void (NAG_CALL *pdeiv)(const Integer *npts, const Integer *npde, const double *t, const double x[], const double y[], double u[]),void (NAG_CALL *monitr)(const Integer *npde, const double *t, const double *dt, const double *dtnew, const logical *tlast, const Integer *nlev, const Integer ngpts[], const double xpts[], const double ypts[], const Integer lsol[], const double sol[], Integer *ierr),const Integer opti[], const double optr[], double rwk[], const Integer *lenrwk, Integer iwk[], const Integer *leniwk, logical lwk[], const Integer *lenlwk, const Integer *itrace, Integer *ind, Integer *ifail)

## 3Description

d03raf integrates the system of PDEs:
 $Fjt,x,y,u,ut,ux,uy,uxx,uxy,uyy=0, j=1,2,…,npde,$ (1)
for $x$ and $y$ in the rectangular domain ${x}_{\mathrm{min}}\le x\le {x}_{\mathrm{max}}$, ${y}_{\mathrm{min}}\le y\le {y}_{\mathrm{max}}$, and time interval ${t}_{0}\le t\le {t}_{\mathrm{out}}$, where the vector $u$ is the set of solution values
 $ux,y,t=u1x,y,t,…,unpdex,y,tT,$
and ${u}_{t}$ denotes partial differentiation with respect to $t$, and similarly for ${u}_{x}$ etc.
The functions ${F}_{j}$ must be supplied by you in pdedef. Similarly the initial values of the functions $u\left(x,y,t\right)$ must be specified at $t={t}_{0}$ in pdeiv.
Note that whilst complete generality is offered by the master equations (1), d03raf is not appropriate for all PDEs. In particular, hyperbolic systems should not be solved using this routine. Also, at least one component of ${u}_{t}$ must appear in the system of PDEs.
The boundary conditions must be supplied by you in bndary in the form
 $Gj t,x,y,u,ut,ux,uy = 0 ,$ (2)
for all $y$ when ${x}_{\mathrm{min}}$ or ${x}_{\mathrm{max}}$ and for all $x$ when $y={y}_{\mathrm{min}}$ or $y={y}_{\mathrm{max}}$ and $j=1,2,\dots ,{\mathbf{npde}}$
The domain is covered by a uniform coarse base grid of size ${n}_{x}×{n}_{y}$ specified by you, and nested finer uniform subgrids are subsequently created in regions with high spatial activity. The refinement is controlled using a space monitor which is computed from the current solution and a user-supplied space tolerance tols. A number of optional parameters, e.g., the maximum number of grid levels at any time, and some weighting factors, can be specified in the arrays opti and optr. Further details of the refinement strategy can be found in Section 9.
The system of PDEs and the boundary conditions are discretized in space on each grid using a standard second-order finite difference scheme (centred on the internal domain and one-sided at the boundaries), and the resulting system of ODEs is integrated in time using a second-order, two-step, implicit BDF method with variable step size. The time integration is controlled using a time monitor computed at each grid level from the current solution and a user-supplied time tolerance tolt, and some further optional user-specified weighting factors held in optr (see Section 9 for details). The time monitor is used to compute a new step size, subject to restrictions on the size of the change between steps, and (optional) user-specified maximum and minimum step sizes held in dt. The step size is adjusted so that the remaining integration interval is an integer number times $\Delta t$. In this way a solution is obtained at $t={t}_{\mathrm{out}}$.
A modified Newton method is used to solve the nonlinear equations arising from the time integration. You may specify (in opti) the maximum number of Newton iterations to be attempted. A Jacobian matrix is calculated at the beginning of each time step. If the Newton process diverges or the maximum number of iterations is exceeded, a new Jacobian is calculated using the most recent iterates and the Newton process is restarted. If convergence is not achieved after the (optional) user-specified maximum number of new Jacobian evaluations, the time step is retried with $\Delta t=\Delta t/4$. The linear systems arising from the Newton iteration are solved using a Bi-CGSTAB iterative method, in combination with ILU preconditioning. The maximum number of iterations can be specified by you in opti.
The solution at all grid levels is stored in the workspace arrays, along with other information needed for a restart (i.e., a continuation call). It is not intended that you extract the solution from these arrays, indeed the necessary information regarding these arrays is not included. The user-supplied monitor monitr should be used to obtain the solution at particular levels and times. monitr is called at the end of every time step, with the last step being identified via the input argument tlast.
Within pdeiv, pdedef, bndary and monitr the data structure is as follows. Each point on a particular grid is given an index (ranging from $1$ to the total number of points on the grid) and all coordinate or solution information is stored in arrays according to this index, e.g., ${\mathbf{x}}\left(i\right)$ and ${\mathbf{y}}\left(i\right)$ contain the $x$- and $y$ coordinate of point $i$, and ${\mathbf{u}}\left(i,j\right)$ contains the $j$th solution component ${u}_{j}$ at point $i$.
Further details of the underlying algorithm can be found in Section 9 and in Blom and Verwer (1993) and Blom et al. (1996) and the references therein.

## 4References

Adjerid S and Flaherty J E (1988) A local refinement finite element method for two-dimensional parabolic systems SIAM J. Sci. Statist. Comput. 9 792–811
Blom J G, Trompert R A and Verwer J G (1996) Algorithm 758. VLUGR2: A vectorizable adaptive grid solver for PDEs in 2D Trans. Math. Software 22 302–328
Blom J G and Verwer J G (1993) VLUGR2: A vectorized local uniform grid refinement code for PDEs in 2D Report NM-R9306 CWI, Amsterdam
Brown P N, Hindmarsh A C and Petzold L R (1994) Using Krylov methods in the solution of large scale differential-algebraic systems SIAM J. Sci. Statist. Comput. 15 1467–1488
Trompert R A (1993) Local uniform grid refinement and systems of coupled partial differential equations Appl. Numer. Maths 12 331–355
Trompert R A and Verwer J G (1993) Analysis of the implicit Euler local uniform grid refinement method SIAM J. Sci. Comput. 14 259–278

## 5Arguments

1:     $\mathbf{npde}$ – IntegerInput
On entry: the number of PDEs in the system.
Constraint: ${\mathbf{npde}}\ge 1$.
2:     $\mathbf{ts}$ – Real (Kind=nag_wp)Input/Output
On entry: the initial value of the independent variable $t$.
On exit: the value of $t$ which has been reached. Normally ${\mathbf{ts}}={\mathbf{tout}}$.
Constraint: ${\mathbf{ts}}<{\mathbf{tout}}$.
3:     $\mathbf{tout}$ – Real (Kind=nag_wp)Input
On entry: the final value of $t$ to which the integration is to be carried out.
4:     $\mathbf{dt}\left(3\right)$ – Real (Kind=nag_wp) arrayInput/Output
On entry: the initial, minimum and maximum time step sizes respectively.
${\mathbf{dt}}\left(1\right)$
Specifies the initial time step size to be used on the first entry, i.e., when ${\mathbf{ind}}=0$. If ${\mathbf{dt}}\left(1\right)=0.0$ then the default value ${\mathbf{dt}}\left(1\right)=0.01×\left({\mathbf{tout}}-{\mathbf{ts}}\right)$ is used. On subsequent entries (${\mathbf{ind}}=1$), the value of ${\mathbf{dt}}\left(1\right)$ is not referenced.
${\mathbf{dt}}\left(2\right)$
Specifies the minimum time step size to be attempted by the integrator. If ${\mathbf{dt}}\left(2\right)=0.0$ the default value  is used.
${\mathbf{dt}}\left(3\right)$
Specifies the maximum time step size to be attempted by the integrator. If ${\mathbf{dt}}\left(3\right)=0.0$ the default value ${\mathbf{dt}}\left(3\right)={\mathbf{tout}}-{\mathbf{ts}}$ is used.
On exit: ${\mathbf{dt}}\left(1\right)$ contains the time step size for the next time step. ${\mathbf{dt}}\left(2\right)$ and ${\mathbf{dt}}\left(3\right)$ are unchanged or set to their default values if zero on entry.
Constraints:
• if ${\mathbf{ind}}=0$, ${\mathbf{dt}}\left(1\right)\ge 0.0$;
• if ${\mathbf{ind}}=0$ and ${\mathbf{dt}}\left(1\right)>0.0$,  and ${\mathbf{dt}}\left(2\right)\le {\mathbf{dt}}\left(1\right)\le {\mathbf{dt}}\left(3\right)$, where the values of ${\mathbf{dt}}\left(2\right)$ and ${\mathbf{dt}}\left(3\right)$ will have been reset to their default values if zero on entry;
• $0\le {\mathbf{dt}}\left(2\right)\le {\mathbf{dt}}\left(3\right)$.
5:     $\mathbf{xmin}$ – Real (Kind=nag_wp)Input
6:     $\mathbf{xmax}$ – Real (Kind=nag_wp)Input
On entry: the extents of the rectangular domain in the $x$-direction, i.e., the $x$ coordinates of the left and right boundaries respectively.
Constraint: ${\mathbf{xmin}}<{\mathbf{xmax}}$ and xmax must be sufficiently distinguishable from xmin for the precision of the machine being used.
7:     $\mathbf{ymin}$ – Real (Kind=nag_wp)Input
8:     $\mathbf{ymax}$ – Real (Kind=nag_wp)Input
On entry: the extents of the rectangular domain in the $y$-direction, i.e., the $y$ coordinates of the lower and upper boundaries respectively.
Constraint: ${\mathbf{ymin}}<{\mathbf{ymax}}$ and ymax must be sufficiently distinguishable from ymin for the precision of the machine being used.
9:     $\mathbf{nx}$ – IntegerInput
On entry: the number of grid points in the $x$-direction (including the boundary points).
Constraint: ${\mathbf{nx}}\ge 4$.
10:   $\mathbf{ny}$ – IntegerInput
On entry: the number of grid points in the $y$-direction (including the boundary points).
Constraint: ${\mathbf{ny}}\ge 4$.
11:   $\mathbf{tols}$ – Real (Kind=nag_wp)Input
On entry: the space tolerance used in the grid refinement strategy ($\sigma$ in equation (4)). See Section 9.2.
Constraint: ${\mathbf{tols}}>0.0$.
12:   $\mathbf{tolt}$ – Real (Kind=nag_wp)Input
On entry: the time tolerance used to determine the time step size ($\tau$ in equation (7)). See Section 9.3.
Constraint: ${\mathbf{tolt}}>0.0$.
13:   $\mathbf{pdedef}$ – Subroutine, supplied by the user.External Procedure
pdedef must evaluate the functions ${F}_{j}$, for $\mathit{j}=1,2,\dots ,{\mathbf{npde}}$, in equation (1) which define the system of PDEs (i.e., the residuals of the resulting ODE system) at all interior points of the domain. Values at points on the boundaries of the domain are ignored and will be overwritten by bndary. pdedef is called for each subgrid in turn.
The specification of pdedef is:
Fortran Interface
 Subroutine pdedef ( npts, npde, t, x, y, u, ut, ux, uy, uxx, uxy, uyy, res)
 Integer, Intent (In) :: npts, npde Real (Kind=nag_wp), Intent (In) :: t, x(npts), y(npts), u(npts,npde), ut(npts,npde), ux(npts,npde), uy(npts,npde), uxx(npts,npde), uxy(npts,npde), uyy(npts,npde) Real (Kind=nag_wp), Intent (Out) :: res(npts,npde)
#include nagmk26.h
 void pdedef (const Integer *npts, const Integer *npde, const double *t, const double x[], const double y[], const double u[], const double ut[], const double ux[], const double uy[], const double uxx[], const double uxy[], const double uyy[], double res[])
1:     $\mathbf{npts}$ – IntegerInput
On entry: the number of grid points in the current grid.
2:     $\mathbf{npde}$ – IntegerInput
On entry: the number of PDEs in the system.
3:     $\mathbf{t}$ – Real (Kind=nag_wp)Input
On entry: the current value of the independent variable $t$.
4:     $\mathbf{x}\left({\mathbf{npts}}\right)$ – Real (Kind=nag_wp) arrayInput
On entry: ${\mathbf{x}}\left(\mathit{i}\right)$ contains the $x$ coordinate of the $\mathit{i}$th grid point, for $\mathit{i}=1,2,\dots ,{\mathbf{npts}}$.
5:     $\mathbf{y}\left({\mathbf{npts}}\right)$ – Real (Kind=nag_wp) arrayInput
On entry: ${\mathbf{y}}\left(\mathit{i}\right)$ contains the $y$ coordinate of the $\mathit{i}$th grid point, for $\mathit{i}=1,2,\dots ,{\mathbf{npts}}$.
6:     $\mathbf{u}\left({\mathbf{npts}},{\mathbf{npde}}\right)$ – Real (Kind=nag_wp) arrayInput
On entry: ${\mathbf{u}}\left(\mathit{i},\mathit{j}\right)$ contains the value of the $\mathit{j}$th PDE component at the $\mathit{i}$th grid point, for $\mathit{i}=1,2,\dots ,{\mathbf{npts}}$ and $\mathit{j}=1,2,\dots ,{\mathbf{npde}}$.
7:     $\mathbf{ut}\left({\mathbf{npts}},{\mathbf{npde}}\right)$ – Real (Kind=nag_wp) arrayInput
On entry: ${\mathbf{ut}}\left(\mathit{i},\mathit{j}\right)$ contains the value of $\frac{\partial u}{\partial t}$ for the $\mathit{j}$th PDE component at the $\mathit{i}$th grid point, for $\mathit{i}=1,2,\dots ,{\mathbf{npts}}$ and $\mathit{j}=1,2,\dots ,{\mathbf{npde}}$.
8:     $\mathbf{ux}\left({\mathbf{npts}},{\mathbf{npde}}\right)$ – Real (Kind=nag_wp) arrayInput
On entry: ${\mathbf{ux}}\left(\mathit{i},\mathit{j}\right)$ contains the value of $\frac{\partial u}{\partial x}$ for the $\mathit{j}$th PDE component at the $\mathit{i}$th grid point, for $\mathit{i}=1,2,\dots ,{\mathbf{npts}}$ and $\mathit{j}=1,2,\dots ,{\mathbf{npde}}$.
9:     $\mathbf{uy}\left({\mathbf{npts}},{\mathbf{npde}}\right)$ – Real (Kind=nag_wp) arrayInput
On entry: ${\mathbf{uy}}\left(\mathit{i},\mathit{j}\right)$ contains the value of $\frac{\partial u}{\partial y}$ for the $\mathit{j}$th PDE component at the $\mathit{i}$th grid point, for $\mathit{i}=1,2,\dots ,{\mathbf{npts}}$ and $\mathit{j}=1,2,\dots ,{\mathbf{npde}}$.
10:   $\mathbf{uxx}\left({\mathbf{npts}},{\mathbf{npde}}\right)$ – Real (Kind=nag_wp) arrayInput
On entry: ${\mathbf{uxx}}\left(\mathit{i},\mathit{j}\right)$ contains the value of $\frac{{\partial }^{2}u}{\partial {x}^{2}}$ for the $\mathit{j}$th PDE component at the $\mathit{i}$th grid point, for $\mathit{i}=1,2,\dots ,{\mathbf{npts}}$ and $\mathit{j}=1,2,\dots ,{\mathbf{npde}}$.
11:   $\mathbf{uxy}\left({\mathbf{npts}},{\mathbf{npde}}\right)$ – Real (Kind=nag_wp) arrayInput
On entry: ${\mathbf{uxy}}\left(\mathit{i},\mathit{j}\right)$ contains the value of $\frac{{\partial }^{2}u}{\partial x\partial y}$ for the $\mathit{j}$th PDE component at the $\mathit{i}$th grid point, for $\mathit{i}=1,2,\dots ,{\mathbf{npts}}$ and $\mathit{j}=1,2,\dots ,{\mathbf{npde}}$.
12:   $\mathbf{uyy}\left({\mathbf{npts}},{\mathbf{npde}}\right)$ – Real (Kind=nag_wp) arrayInput
On entry: ${\mathbf{uyy}}\left(\mathit{i},\mathit{j}\right)$ contains the value of $\frac{{\partial }^{2}u}{\partial {y}^{2}}$ for the $\mathit{j}$th PDE component at the $\mathit{i}$th grid point, for $\mathit{i}=1,2,\dots ,{\mathbf{npts}}$ and $\mathit{j}=1,2,\dots ,{\mathbf{npde}}$.
13:   $\mathbf{res}\left({\mathbf{npts}},{\mathbf{npde}}\right)$ – Real (Kind=nag_wp) arrayOutput
On exit: ${\mathbf{res}}\left(\mathit{i},\mathit{j}\right)$ must contain the value of ${F}_{\mathit{j}}$, for $\mathit{j}=1,2,\dots ,{\mathbf{npde}}$, at the $\mathit{i}$th grid point, for $\mathit{i}=1,2,\dots ,{\mathbf{npts}}$, although the residuals at boundary points will be ignored (and overwritten later on) and so they need not be specified here.
pdedef must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which d03raf is called. Arguments denoted as Input must not be changed by this procedure.
Note: pdedef should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by d03raf. If your code inadvertently does return any NaNs or infinities, d03raf is likely to produce unexpected results.
14:   $\mathbf{bndary}$ – Subroutine, supplied by the user.External Procedure
bndary must evaluate the functions ${G}_{\mathit{j}}$, for $\mathit{j}=1,2,\dots ,{\mathbf{npde}}$, in equation (2) which define the boundary conditions at all boundary points of the domain. Residuals at interior points must not be altered by this subroutine.
The specification of bndary is:
Fortran Interface
 Subroutine bndary ( npts, npde, t, x, y, u, ut, ux, uy, lbnd, res)
 Integer, Intent (In) :: npts, npde, nbpts, lbnd(nbpts) Real (Kind=nag_wp), Intent (In) :: t, x(npts), y(npts), u(npts,npde), ut(npts,npde), ux(npts,npde), uy(npts,npde) Real (Kind=nag_wp), Intent (Inout) :: res(npts,npde)
#include nagmk26.h
 void bndary (const Integer *npts, const Integer *npde, const double *t, const double x[], const double y[], const double u[], const double ut[], const double ux[], const double uy[], const Integer *nbpts, const Integer lbnd[], double res[])
1:     $\mathbf{npts}$ – IntegerInput
On entry: the number of grid points in the current grid.
2:     $\mathbf{npde}$ – IntegerInput
On entry: the number of PDEs in the system.
3:     $\mathbf{t}$ – Real (Kind=nag_wp)Input
On entry: the current value of the independent variable $t$.
4:     $\mathbf{x}\left({\mathbf{npts}}\right)$ – Real (Kind=nag_wp) arrayInput
On entry: ${\mathbf{x}}\left(\mathit{i}\right)$ contains the $x$ coordinate of the $\mathit{i}$th grid point, for $\mathit{i}=1,2,\dots ,{\mathbf{npts}}$.
5:     $\mathbf{y}\left({\mathbf{npts}}\right)$ – Real (Kind=nag_wp) arrayInput
On entry: ${\mathbf{y}}\left(\mathit{i}\right)$ contains the $y$ coordinate of the $\mathit{i}$th grid point, for $\mathit{i}=1,2,\dots ,{\mathbf{npts}}$.
6:     $\mathbf{u}\left({\mathbf{npts}},{\mathbf{npde}}\right)$ – Real (Kind=nag_wp) arrayInput
On entry: ${\mathbf{u}}\left(\mathit{i},\mathit{j}\right)$ contains the value of the $\mathit{j}$th PDE component at the $\mathit{i}$th grid point, for $\mathit{i}=1,2,\dots ,{\mathbf{npts}}$ and $\mathit{j}=1,2,\dots ,{\mathbf{npde}}$.
7:     $\mathbf{ut}\left({\mathbf{npts}},{\mathbf{npde}}\right)$ – Real (Kind=nag_wp) arrayInput
On entry: ${\mathbf{ut}}\left(\mathit{i},\mathit{j}\right)$ contains the value of $\frac{\partial u}{\partial t}$ for the $\mathit{j}$th PDE component at the $\mathit{i}$th grid point, for $\mathit{i}=1,2,\dots ,{\mathbf{npts}}$ and $\mathit{j}=1,2,\dots ,{\mathbf{npde}}$.
8:     $\mathbf{ux}\left({\mathbf{npts}},{\mathbf{npde}}\right)$ – Real (Kind=nag_wp) arrayInput
On entry: ${\mathbf{ux}}\left(\mathit{i},\mathit{j}\right)$ contains the value of $\frac{\partial u}{\partial x}$ for the $\mathit{j}$th PDE component at the $\mathit{i}$th grid point, for $\mathit{i}=1,2,\dots ,{\mathbf{npts}}$ and $\mathit{j}=1,2,\dots ,{\mathbf{npde}}$.
9:     $\mathbf{uy}\left({\mathbf{npts}},{\mathbf{npde}}\right)$ – Real (Kind=nag_wp) arrayInput
On entry: ${\mathbf{uy}}\left(\mathit{i},\mathit{j}\right)$ contains the value of $\frac{\partial u}{\partial y}$ for the $\mathit{j}$th PDE component at the $\mathit{i}$th grid point, for $\mathit{i}=1,2,\dots ,{\mathbf{npts}}$ and $\mathit{j}=1,2,\dots ,{\mathbf{npde}}$.
10:   $\mathbf{nbpts}$ – IntegerInput
On entry: the number of boundary points in the grid.
11:   $\mathbf{lbnd}\left({\mathbf{nbpts}}\right)$ – Integer arrayInput
On entry: ${\mathbf{lbnd}}\left(\mathit{i}\right)$ contains the grid index for the $\mathit{i}$th boundary point, for $\mathit{i}=1,2,\dots ,{\mathbf{nbpts}}$. Hence the $\mathit{i}$th boundary point has coordinates ${\mathbf{x}}\left({\mathbf{lbnd}}\left(\mathit{i}\right)\right)$ and ${\mathbf{y}}\left({\mathbf{lbnd}}\left(\mathit{i}\right)\right)$, and the corresponding solution values are ${\mathbf{u}}\left({\mathbf{lbnd}}\left(\mathit{i}\right),{\mathbf{npde}}\right)$, etc.
12:   $\mathbf{res}\left({\mathbf{npts}},{\mathbf{npde}}\right)$ – Real (Kind=nag_wp) arrayInput/Output
On entry: ${\mathbf{res}}\left(\mathit{i},\mathit{j}\right)$ contains the value of ${F}_{\mathit{j}}$, for $\mathit{i}=1,2,\dots ,{\mathbf{npde}}$, at the $\mathit{i}$th grid point, for $\mathit{i}=1,2,\dots ,{\mathbf{npts}}$, as returned by pdedef. The residuals at the boundary points will be overwritten and so need not have been set by pdedef.
On exit: ${\mathbf{res}}\left({\mathbf{lbnd}}\left(\mathit{i}\right),\mathit{j}\right)$ must contain the value of ${G}_{\mathit{j}}$, for $\mathit{j}=1,2,\dots ,{\mathbf{npde}}$, at the $\mathit{i}$th boundary point, for $\mathit{i}=1,2,\dots ,{\mathbf{nbpts}}$.
Note: elements of res corresponding to interior points must not be altered.
bndary must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which d03raf is called. Arguments denoted as Input must not be changed by this procedure.
Note: bndary should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by d03raf. If your code inadvertently does return any NaNs or infinities, d03raf is likely to produce unexpected results.
15:   $\mathbf{pdeiv}$ – Subroutine, supplied by the user.External Procedure
pdeiv must specify the initial values of the PDE components $u$ at all points in the grid. pdeiv is not referenced if, on entry, ${\mathbf{ind}}=1$.
The specification of pdeiv is:
Fortran Interface
 Subroutine pdeiv ( npts, npde, t, x, y, u)
 Integer, Intent (In) :: npts, npde Real (Kind=nag_wp), Intent (In) :: t, x(npts), y(npts) Real (Kind=nag_wp), Intent (Out) :: u(npts,npde)
#include nagmk26.h
 void pdeiv (const Integer *npts, const Integer *npde, const double *t, const double x[], const double y[], double u[])
1:     $\mathbf{npts}$ – IntegerInput
On entry: the number of grid points in the grid.
2:     $\mathbf{npde}$ – IntegerInput
On entry: the number of PDEs in the system.
3:     $\mathbf{t}$ – Real (Kind=nag_wp)Input
On entry: the (initial) value of the independent variable $t$.
4:     $\mathbf{x}\left({\mathbf{npts}}\right)$ – Real (Kind=nag_wp) arrayInput
On entry: ${\mathbf{x}}\left(\mathit{i}\right)$ contains the $x$ coordinate of the $\mathit{i}$th grid point, for $\mathit{i}=1,2,\dots ,{\mathbf{npts}}$.
5:     $\mathbf{y}\left({\mathbf{npts}}\right)$ – Real (Kind=nag_wp) arrayInput
On entry: ${\mathbf{y}}\left(\mathit{i}\right)$ contains the $y$ coordinate of the $\mathit{i}$th grid point, for $\mathit{i}=1,2,\dots ,{\mathbf{npts}}$.
6:     $\mathbf{u}\left({\mathbf{npts}},{\mathbf{npde}}\right)$ – Real (Kind=nag_wp) arrayOutput
On exit: ${\mathbf{u}}\left(\mathit{i},\mathit{j}\right)$ must contain the value of the $\mathit{j}$th PDE component at the $\mathit{i}$th grid point, for $\mathit{i}=1,2,\dots ,{\mathbf{npts}}$ and $\mathit{j}=1,2,\dots ,{\mathbf{npde}}$.
pdeiv must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which d03raf is called. Arguments denoted as Input must not be changed by this procedure.
Note: pdeiv should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by d03raf. If your code inadvertently does return any NaNs or infinities, d03raf is likely to produce unexpected results.
16:   $\mathbf{monitr}$ – Subroutine, supplied by the user.External Procedure
monitr is called by d03raf at the end of every successful time step, and may be used to examine or print the solution or perform other tasks such as error calculations, particularly at the final time step, indicated by the argument tlast. The input arguments contain information about the grid and solution at all grid levels used.
monitr can also be used to force an immediate tidy termination of the solution process and return to the calling program.
The specification of monitr is:
Fortran Interface
 Subroutine monitr ( npde, t, dt, nlev, xpts, ypts, lsol, sol, ierr)
 Integer, Intent (In) :: npde, nlev, ngpts(nlev), lsol(nlev) Integer, Intent (Inout) :: ierr Real (Kind=nag_wp), Intent (In) :: t, dt, dtnew, xpts(*), ypts(*), sol(*) Logical, Intent (In) :: tlast
#include nagmk26.h
 void monitr (const Integer *npde, const double *t, const double *dt, const double *dtnew, const logical *tlast, const Integer *nlev, const Integer ngpts[], const double xpts[], const double ypts[], const Integer lsol[], const double sol[], Integer *ierr)
1:     $\mathbf{npde}$ – IntegerInput
On entry: the number of PDEs in the system.
2:     $\mathbf{t}$ – Real (Kind=nag_wp)Input
On entry: the current value of the independent variable $t$, i.e., the time at the end of the integration step just completed.
3:     $\mathbf{dt}$ – Real (Kind=nag_wp)Input
On entry: the current time step size $\Delta t$, i.e., the time step size used for the integration step just completed.
4:     $\mathbf{dtnew}$ – Real (Kind=nag_wp)Input
On entry: the step size that will be used for the next time step.
5:     $\mathbf{tlast}$ – LogicalInput
On entry: indicates if intermediate or final time step. ${\mathbf{tlast}}=\mathrm{.FALSE.}$ for an intermediate step, ${\mathbf{tlast}}=\mathrm{.TRUE.}$ for the last call to monitr before returning to your program.
6:     $\mathbf{nlev}$ – IntegerInput
On entry: the number of grid levels used at time t.
7:     $\mathbf{ngpts}\left({\mathbf{nlev}}\right)$ – Integer arrayInput
On entry: ${\mathbf{ngpts}}\left(\mathit{l}\right)$ contains the number of grid points at level $\mathit{l}$, for $\mathit{l}=1,2,\dots ,{\mathbf{nlev}}$.
8:     $\mathbf{xpts}\left(*\right)$ – Real (Kind=nag_wp) arrayInput
On entry: contains the $x$ coordinates of the grid points in each level in turn, i.e., ${\mathbf{x}}\left(\mathit{i}\right)$, for $\mathit{i}=1,2,\dots ,{\mathbf{ngpts}}\left(\mathit{l}\right)$ and $\mathit{l}=1,2,\dots ,{\mathbf{nlev}}$.
So for level $\mathit{l}$, ${\mathbf{x}}\left(\mathit{i}\right)={\mathbf{xpts}}\left(\mathit{k}+\mathit{i}\right)$, where $\mathit{k}={\mathbf{ngpts}}\left(1\right)+{\mathbf{ngpts}}\left(2\right)+\cdots +{\mathbf{ngpts}}\left(\mathit{l}-1\right)$, for $\mathit{i}=1,2,\dots ,{\mathbf{ngpts}}\left(\mathit{l}\right)$ and $\mathit{l}=1,2,\dots ,{\mathbf{nlev}}$.
9:     $\mathbf{ypts}\left(*\right)$ – Real (Kind=nag_wp) arrayInput
On entry: contains the $y$ coordinates of the grid points in each level in turn, i.e., ${\mathbf{y}}\left(\mathit{i}\right)$, for $\mathit{i}=1,2,\dots ,{\mathbf{ngpts}}\left(\mathit{l}\right)$ and $\mathit{l}=1,2,\dots ,{\mathbf{nlev}}$.
So for level $\mathit{l}$, ${\mathbf{y}}\left(\mathit{i}\right)={\mathbf{ypts}}\left(\mathit{k}+\mathit{i}\right)$, where $\mathit{k}={\mathbf{ngpts}}\left(1\right)+{\mathbf{ngpts}}\left(2\right)+\cdots +{\mathbf{ngpts}}\left(\mathit{l}-1\right)$, for $\mathit{i}=1,2,\dots ,{\mathbf{ngpts}}\left(\mathit{l}\right)$ and $\mathit{l}=1,2,\dots ,{\mathbf{nlev}}$.
10:   $\mathbf{lsol}\left({\mathbf{nlev}}\right)$ – Integer arrayInput
On entry: ${\mathbf{lsol}}\left(\mathit{l}\right)$ contains the pointer to the solution in sol at grid level $\mathit{l}$ and time t. (${\mathbf{lsol}}\left(\mathit{l}\right)$ actually contains the array index immediately preceding the start of the solution in sol.)
11:   $\mathbf{sol}\left(*\right)$ – Real (Kind=nag_wp) arrayInput
On entry: contains the solution ${\mathbf{u}}\left({\mathbf{ngpts}}\left(\mathit{l}\right),{\mathbf{npde}}\right)$ at time t for each grid level $\mathit{l}$ in turn, positioned according to lsol, i.e., for level $\mathit{l}$, ${\mathbf{u}}\left(\mathit{i},\mathit{j}\right)={\mathbf{sol}}\left({\mathbf{lsol}}\left(\mathit{l}\right)+\left(\mathit{j}-1\right)×{\mathbf{ngpts}}\left(\mathit{l}\right)+\mathit{i}\right)$, for $\mathit{i}=1,2,\dots ,{\mathbf{ngpts}}\left(\mathit{l}\right)$, $\mathit{j}=1,2,\dots ,{\mathbf{npde}}$ and $\mathit{l}=1,2,\dots ,{\mathbf{nlev}}$.
12:   $\mathbf{ierr}$ – IntegerInput/Output
On entry: will be set to $0$.
On exit: should be set to $1$ to force a tidy termination and an immediate return to the calling program with ${\mathbf{ifail}}={\mathbf{4}}$. ierr should remain unchanged otherwise.
monitr must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which d03raf is called. Arguments denoted as Input must not be changed by this procedure.
17:   $\mathbf{opti}\left(4\right)$ – Integer arrayInput
On entry: may be set to control various options available in the integrator.
${\mathbf{opti}}\left(1\right)=0$
All the default options are employed.
${\mathbf{opti}}\left(1\right)>0$
The default value of ${\mathbf{opti}}\left(\mathit{i}\right)$, for $\mathit{i}=2,3,4$, can be obtained by setting ${\mathbf{opti}}\left(\mathit{i}\right)=0$.
${\mathbf{opti}}\left(1\right)$
Specifies the maximum number of grid levels allowed (including the base grid). ${\mathbf{opti}}\left(1\right)\ge 0$. The default value is ${\mathbf{opti}}\left(1\right)=3$.
${\mathbf{opti}}\left(2\right)$
Specifies the maximum number of Jacobian evaluations allowed during each nonlinear equations solution. ${\mathbf{opti}}\left(2\right)\ge 0$. The default value is ${\mathbf{opti}}\left(2\right)=2$.
${\mathbf{opti}}\left(3\right)$
Specifies the maximum number of Newton iterations in each nonlinear equations solution. ${\mathbf{opti}}\left(3\right)\ge 0$. The default value is ${\mathbf{opti}}\left(3\right)=10$.
${\mathbf{opti}}\left(4\right)$
Specifies the maximum number of iterations in each linear equations solution. ${\mathbf{opti}}\left(4\right)\ge 0$. The default value is ${\mathbf{opti}}\left(4\right)=100$.
Constraint: ${\mathbf{opti}}\left(1\right)\ge 0$ and if ${\mathbf{opti}}\left(1\right)>0$, ${\mathbf{opti}}\left(\mathit{i}\right)\ge 0$, for $\mathit{i}=2,3,4$.
18:   $\mathbf{optr}\left(3,{\mathbf{npde}}\right)$ – Real (Kind=nag_wp) arrayInput
On entry: may be used to specify the optional vectors ${u}^{\mathrm{max}}$, ${w}^{\mathrm{s}}$ and ${w}^{\mathrm{t}}$ in the space and time monitors (see Section 9).
If an optional vector is not required then all its components should be set to $1.0$.
${\mathbf{optr}}\left(1,\mathit{j}\right)$, for $\mathit{j}=1,2,\dots ,{\mathbf{npde}}$, specifies ${u}_{\mathit{j}}^{\mathrm{max}}$, the approximate maximum absolute value of the $\mathit{j}$th component of $u$, as used in (4) and (7). ${\mathbf{optr}}\left(1,\mathit{j}\right)>0.0$, for $\mathit{j}=1,2,\dots ,{\mathbf{npde}}$.
${\mathbf{optr}}\left(2,\mathit{j}\right)$, for $\mathit{j}=1,2,\dots ,{\mathbf{npde}}$, specifies ${w}_{\mathit{j}}^{\mathrm{s}}$, the weighting factors used in the space monitor (see (4)) to indicate the relative importance of the $\mathit{j}$th component of $u$ on the space monitor. ${\mathbf{optr}}\left(2,\mathit{j}\right)\ge 0.0$, for $\mathit{j}=1,2,\dots ,{\mathbf{npde}}$.
${\mathbf{optr}}\left(3,\mathit{j}\right)$, for $\mathit{j}=1,2,\dots ,{\mathbf{npde}}$, specifies ${w}_{\mathit{j}}^{\mathrm{t}}$, the weighting factors used in the time monitor (see (6)) to indicate the relative importance of the $\mathit{j}$th component of $u$ on the time monitor. ${\mathbf{optr}}\left(3,\mathit{j}\right)\ge 0.0$, for $\mathit{j}=1,2,\dots ,{\mathbf{npde}}$.
Constraints:
• ${\mathbf{optr}}\left(1,\mathit{j}\right)>0.0$, for $\mathit{j}=1,2,\dots ,{\mathbf{npde}}$;
• ${\mathbf{optr}}\left(\mathit{i},\mathit{j}\right)\ge 0.0$, for $\mathit{i}=2,3$ and $\mathit{j}=1,2,\dots ,{\mathbf{npde}}$.
19:   $\mathbf{rwk}\left({\mathbf{lenrwk}}\right)$ – Real (Kind=nag_wp) arrayCommunication Array
20:   $\mathbf{lenrwk}$ – IntegerInput
On entry: the dimension of the array rwk as declared in the (sub)program from which d03raf is called.
The required value of lenrwk cannot be determined exactly in advance, but a suggested value is
 $lenrwk=maxpts×npde×5×l+18×npde+9+2×maxpts,$
where $\mathit{l}={\mathbf{opti}}\left(1\right)$ if ${\mathbf{opti}}\left(1\right)\ne 0$ and $\mathit{l}=3$ otherwise, and $\mathit{maxpts}$ is the expected maximum number of grid points at any one level. If during the execution the supplied value is found to be too small then the routine returns with ${\mathbf{ifail}}={\mathbf{3}}$ and an estimated required size is printed on the current error message unit (see x04aaf).
Constraint: ${\mathbf{lenrwk}}\ge {\mathbf{nx}}×{\mathbf{ny}}×{\mathbf{npde}}×\left(14+18×{\mathbf{npde}}\right)+2×{\mathbf{nx}}×{\mathbf{ny}}$ (the required size for the initial grid).
21:   $\mathbf{iwk}\left({\mathbf{leniwk}}\right)$ – Integer arrayCommunication Array
On entry: if ${\mathbf{ind}}=0$, iwk need not be set. Otherwise iwk must remain unchanged from a previous call to d03raf.
On exit: the following components of the array iwk concern the efficiency of the integration. Here, $m$ is the maximum number of grid levels allowed ($m={\mathbf{opti}}\left(1\right)$ if ${\mathbf{opti}}\left(1\right)>1$ and $m=3$ otherwise), and $\mathit{l}$ is a grid level taking the values $\mathit{l}=1,2,\dots ,\mathit{nl}$, where $\mathit{nl}$ is the number of levels used.
${\mathbf{iwk}}\left(1\right)$
Contains the number of steps taken in time.
${\mathbf{iwk}}\left(2\right)$
Contains the number of rejected time steps.
${\mathbf{iwk}}\left(2+\mathit{l}\right)$
Contains the total number of residual evaluations performed (i.e., the number of times pdedef was called) at grid level $\mathit{l}$.
${\mathbf{iwk}}\left(2+m+\mathit{l}\right)$
Contains the total number of Jacobian evaluations performed at grid level $\mathit{l}$.
${\mathbf{iwk}}\left(2+2×m+\mathit{l}\right)$
Contains the total number of Newton iterations performed at grid level $\mathit{l}$.
${\mathbf{iwk}}\left(2+3×m+\mathit{l}\right)$
Contains the total number of linear solver iterations performed at grid level $\mathit{l}$.
${\mathbf{iwk}}\left(2+4×m+\mathit{l}\right)$
Contains the maximum number of Newton iterations performed at any one time step at grid level $\mathit{l}$.
${\mathbf{iwk}}\left(2+5×m+\mathit{l}\right)$
Contains the maximum number of linear solver iterations performed at any one time step at grid level $\mathit{l}$.
Note: the total and maximum numbers are cumulative over all calls to d03raf. If the specified maximum number of Newton or linear solver iterations is exceeded at any stage, the maximums above are set to the specified maximum plus one.
22:   $\mathbf{leniwk}$ – IntegerInput
On entry: the dimension of the array iwk as declared in the (sub)program from which d03raf is called.
The required value of leniwk cannot be determined exactly in advance, but a suggested value is
 $leniwk=maxpts×14+5×m+7×m+2,$
where $\mathit{maxpts}$ is the expected maximum number of grid points at any one level and $m={\mathbf{opti}}\left(1\right)$ if ${\mathbf{opti}}\left(1\right)>0$ and $m=3$ otherwise. If during the execution the supplied value is found to be too small then the routine returns with ${\mathbf{ifail}}={\mathbf{3}}$ and an estimated required size is printed on the current error message unit (see x04aaf).
Constraint: ${\mathbf{leniwk}}\ge 19×{\mathbf{nx}}×{\mathbf{ny}}+9$ (the required size for the initial grid).
23:   $\mathbf{lwk}\left({\mathbf{lenlwk}}\right)$ – Logical arrayWorkspace
24:   $\mathbf{lenlwk}$ – IntegerInput
On entry: the dimension of the array lwk as declared in the (sub)program from which d03raf is called.
The required value of lenlwk cannot be determined exactly in advanced, but a suggested value is
 $lenlwk=maxpts+1,$
where $\mathit{maxpts}$ is the expected maximum number of grid points at any one level. If during the execution the supplied value is found to be too small then the routine returns with ${\mathbf{ifail}}={\mathbf{3}}$ and an estimated required size is printed on the current error message unit (see x04aaf).
Constraint: ${\mathbf{lenlwk}}\ge {\mathbf{nx}}×{\mathbf{ny}}+1$ (the required size for the initial grid).
25:   $\mathbf{itrace}$ – IntegerInput
On entry: the level of trace information required from d03raf. itrace may take the value $-1$, $0$, $1$, $2$ or $3$.
${\mathbf{itrace}}=-1$
No output is generated.
${\mathbf{itrace}}=0$
Only warning messages are printed.
${\mathbf{itrace}}>0$
Output from the underlying solver is printed on the current advisory message unit (see x04abf). This output contains details of the time integration, the nonlinear iteration and the linear solver.
If ${\mathbf{itrace}}<-1$, $-1$ is assumed and similarly if ${\mathbf{itrace}}>3$, $3$ is assumed.
The advisory messages are given in greater detail as itrace increases. Setting ${\mathbf{itrace}}=1$ allows you to monitor the progress of the integration without possibly excessive information.
26:   $\mathbf{ind}$ – IntegerInput/Output
On entry: must be set to $0$ or $1$, alternatively $10$ or $11$.
${\mathbf{ind}}=0$
Starts the integration in time. pdedef is assumed to be serial.
${\mathbf{ind}}=1$
Continues the integration after an earlier exit from the routine. In this case, only the following parameters may be reset between calls to d03raf: tout, dt, tols, tolt, opti, optr, itrace and ifail. pdedef is assumed to be serial.
${\mathbf{ind}}=10$
Starts the integration in time. pdedef is assumed to have been parallelized by you, as described in Section 8. In all other respects, this is equivalent to ${\mathbf{ind}}=0$.
${\mathbf{ind}}=11$
Continues the integration after an earlier exit from the routine. In this case, only the following parameters may be reset between calls to d03raf: tout, dt, tols, tolt, opti, optr, itrace and ifail. pdedef is assumed to have been parallelized by you, as described in Section 8. In all other respects, this is equivalent to ${\mathbf{ind}}=1$.
Constraint: $0\le {\mathbf{ind}}\le 1$ or $10\le {\mathbf{ind}}\le 11$.
On exit: ${\mathbf{ind}}=1$, if ind on input was $0$ or $1$, or ${\mathbf{ind}}=11$, if ind on input was $10$ or $11$.
Note: for users of serial versions of the NAG Library, it is recommended that you only use ${\mathbf{ind}}=0$ or $1$. See Section 8 for more information on the use of ind.
27:   $\mathbf{ifail}$ – IntegerInput/Output
On entry: ifail must be set to $0$, $-1\text{​ or ​}1$. If you are unfamiliar with this argument you should refer to Section 3.4 in How to Use the NAG Library and its Documentation for details.
For environments where it might be inappropriate to halt program execution when an error is detected, the value $-1\text{​ or ​}1$ is recommended. If the output of error messages is undesirable, then the value $1$ is recommended. Otherwise, if you are not familiar with this argument, the recommended value is $0$. When the value $-\mathbf{1}\text{​ or ​}\mathbf{1}$ is used it is essential to test the value of ifail on exit.
On exit: ${\mathbf{ifail}}={\mathbf{0}}$ unless the routine detects an error or a warning has been flagged (see Section 6).

## 6Error Indicators and Warnings

If on entry ${\mathbf{ifail}}=0$ or $-1$, explanatory error messages are output on the current error message unit (as defined by x04aaf).
Errors or warnings detected by the routine:
${\mathbf{ifail}}=1$
 On entry, ${\mathbf{npde}}<1$, or ${\mathbf{tout}}\le {\mathbf{ts}}$, or tout is too close to ts, or ${\mathbf{ind}}=0$ and ${\mathbf{dt}}\left(1\right)<0.0$, or ${\mathbf{dt}}\left(i\right)<0.0$, for $\mathit{i}=2\text{​ or ​}3$, or ${\mathbf{dt}}\left(2\right)>{\mathbf{dt}}\left(3\right)$, or ${\mathbf{ind}}=0$ and , or ${\mathbf{ind}}=0$ and ${\mathbf{dt}}\left(1\right)>{\mathbf{tout}}-{\mathbf{ts}}$, or ${\mathbf{ind}}=0$ and ${\mathbf{dt}}\left(1\right)<{\mathbf{dt}}\left(2\right)$ or ${\mathbf{dt}}\left(1\right)>{\mathbf{dt}}\left(3\right)$, or ${\mathbf{xmin}}\ge {\mathbf{xmax}}$, or xmax too close to xmin, or ${\mathbf{ymin}}\ge {\mathbf{ymax}}$, or ymax too close to ymin, or nx or ${\mathbf{ny}}<4$, or tols or ${\mathbf{tolt}}\le 0.0$, or ${\mathbf{opti}}\left(1\right)<0$, or ${\mathbf{opti}}\left(1\right)>0$ and ${\mathbf{opti}}\left(j\right)<0$, for $\mathit{j}=2$, $3$ or $4$, or ${\mathbf{optr}}\left(1,j\right)\le 0.0$, for some $j=1,2,\dots ,{\mathbf{npde}}$, or ${\mathbf{optr}}\left(2,j\right)<0.0$, for some $j=1,2,\dots ,{\mathbf{npde}}$, or ${\mathbf{optr}}\left(3,j\right)<0.0$, for some $j=1,2,\dots ,{\mathbf{npde}}$, or lenrwk, leniwk or lenlwk too small for initial grid level, or ${\mathbf{ind}}\ne 0$ or $1$, or ${\mathbf{ind}}=1$ on initial entry to d03raf.
${\mathbf{ifail}}=2$
The time step size to be attempted is less than the specified minimum size. This may occur following time step failures and subsequent step size reductions caused by one or more of the following:
• the requested accuracy could not be achieved, i.e., tolt is too small,
• the maximum number of linear solver iterations, Newton iterations or Jacobian evaluations is too small,
• ILU decomposition of the Jacobian matrix could not be performed, possibly due to singularity of the Jacobian.
Setting itrace to a higher value may provide further information.
In the latter two cases you are advised to check their problem formulation in pdedef and/or bndary, and the initial values in pdeiv if appropriate.
${\mathbf{ifail}}=3$
One or more of the workspace arrays is too small for the required number of grid points. An estimate of the required sizes for the current stage is output, but more space may be required at a later stage.
${\mathbf{ifail}}=4$
ierr was set to $1$ in monitr, forcing control to be passed back to calling program. Integration was successful as far as ${\mathbf{t}}={\mathbf{ts}}$.
${\mathbf{ifail}}=5$
The integration has been completed but the maximum number of levels specified in ${\mathbf{opti}}\left(1\right)$ was insufficient at one or more time steps, meaning that the requested space accuracy could not be achieved. To avoid this warning either increase the value of ${\mathbf{opti}}\left(1\right)$ or decrease the value of tols.
${\mathbf{ifail}}=-99$
See Section 3.9 in How to Use the NAG Library and its Documentation for further information.
${\mathbf{ifail}}=-399$
Your licence key may have expired or may not have been installed correctly.
See Section 3.8 in How to Use the NAG Library and its Documentation for further information.
${\mathbf{ifail}}=-999$
Dynamic memory allocation failed.
See Section 3.7 in How to Use the NAG Library and its Documentation for further information.

## 7Accuracy

There are three sources of error in the algorithm: space and time discretization, and interpolation (linear) between grid levels. The space and time discretization errors are controlled separately using the arguments tols and tolt described in the following section, and you should test the effects of varying these arguments. Interpolation errors are generally implicitly controlled by the refinement criterion since in areas where interpolation errors are potentially large, the space monitor will also be large. It can be shown that the global spatial accuracy is comparable to that which would be obtained on a uniform grid of the finest grid size. A full error analysis can be found in Trompert and Verwer (1993).

## 8Parallelism and Performance

d03raf is not thread safe and should not be called from a multithreaded user program. Please see Section 3.12.1 in How to Use the NAG Library and its Documentation for more information on thread safety.
d03raf is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
d03raf makes calls to BLAS and/or LAPACK routines, which may be threaded within the vendor library used by this implementation. Consult the documentation for the vendor library for further information.
Please consult the X06 Chapter Introduction for information on how to control and interrogate the OpenMP environment used within this routine. Please also consult the Users' Note for your implementation for any additional implementation-specific information.
d03raf requires a user-supplied routine pdedef to evaluate the functions ${F}_{j}$, for $\mathit{j}=1,2,\dots ,{\mathbf{npde}}$. The parallelism within d03raf will be more efficient if pdedef can also be parallelized. This is often the case, but you must add some OpenMP directives to your version of pdedef to implement the parallelism. For example, if the body of code for pdedef is as follows (adapted from the first test case in the document for d03raf):
```res(1:npts,1:npde) = ut(1:npts,1:npde) - diffusion*(uxx(1:npts,1:         &
npde)+uyy(1:npts,1:npde)) - damkohler*(one+heat_release-u(1:npts,       &
1:npde))*exp(-activ_energy/u(1:npts,1:npde))```
This example can be parallelized, as the updating of res for each value in the range $1,\dots ,{\mathbf{npts}}$ is independent of every other value. Thus this should be parallelized in OpenMP (using an explicit loop rather than Fortran array syntax) as follows:
```!\$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```
Note that the OpenMP PARALLEL directive must not be specified, as the OpenMP DO directive will bind to the PARALLEL region within the d03raf code. Also note that this assumes the default OpenMP behaviour that all variables are SHARED, except for loop indices that are PRIVATE.
To avoid problems for existing library users, who will not have specified any OpenMP directives in their pdedef routine, the default assumption of d03raf is that pdedef has not been parallelized, and executes calls to pdedef in serial mode. You must indicate that pdedef has been parallelized by setting ind to $10$ or $11$ as appropriate. See Section 5 for details.
If the code within pdedef cannot be parallelized, you must not add any OpenMP directives to your code, and must not set ind to $10$ or $11$. If ind is set to $10$ or $11$ and pdedef has not been parallelized, results on multiple threads will be unpredictable and may give rise to incorrect results and/or program crashes or deadlocks. Please contact NAG for advice if required. Overloading ind in this manner is not entirely satisfactory, consequently it is likely that replacement interfaces for d03raf will be included in a future NAG Library release.

### 9.1Algorithm Outline

The local uniform grid refinement method is summarised as follows:
1. Initialize the course base grid, an initial solution and an initial time step.
2. Solve the system of PDEs on the current grid with the current time step.
3. If the required accuracy in space and the maximum number of grid levels have not yet been reached:
 (a) Determine new finer grid at forward time level. (b) Get solution values at previous time level(s) on new grid. (c) Interpolate internal boundary values from old grid at forward time. (d) Get initial values for the Newton process at forward time. (e) Go to $2$.
4. Update the coarser grid solution using the finer grid values.
5. Estimate error in time integration. If time error is acceptable advance time level.
6. Determine new step size then go to $2$ with coarse base as current grid.

### 9.2Refinement Strategy

For each grid point $i$ a space monitor ${\mu }_{i}^{\mathrm{s}}$ is determined by
 $μ i s = max j=1,npde γ j Δ x 2 ∂ 2 ∂ x 2 u j x i , y i ,t + Δ y 2 ∂ 2 ∂ y 2 u j x i , y i ,t ,$ (3)
where $\Delta x$ and $\Delta y$ are the grid widths in the $x$ and $y$ directions; and ${x}_{i}$, ${y}_{i}$ are the $x$ and $y$ coordinates at grid point $i$. The argument ${\gamma }_{j}$ is obtained from
 $γj=wjs ujmax σ ,$ (4)
where $\sigma$ is the user-supplied space tolerance; ${w}_{j}^{s}$ is a weighting factor for the relative importance of the $j$th PDE component on the space monitor; and ${u}_{j}^{\mathrm{max}}$ is the approximate maximum absolute value of the $j$th component. A value for $\sigma$ must be supplied by you. Values for ${w}_{j}^{\mathrm{s}}$ and ${u}_{j}^{\mathrm{max}}$ must also be supplied but may be set to the value $1.0$ if little information about the solution is known.
A new level of refinement is created if
 $maxiμis>0.9 or 1.0,$ (5)
depending on the grid level at the previous step in order to avoid fluctuations in the number of grid levels between time steps. If (5) is satisfied then all grid points for which ${\mu }_{i}^{s}>0.25$ are flagged and surrounding cells are quartered in size.
No derefinement takes place as such, since at each time step the solution on the base grid is computed first and new finer grids are then created based on the new solution. Hence derefinement occurs implicitly. See Section 9.1.

### 9.3Time Integration

The time integration is controlled using a time monitor calculated at each level $l$ up to the maximum level used, given by
 $μlt=1N∑j=1npdewjt∑i=1 ngptsl Δtαij utxi,yi,t 2$ (6)
where ${\mathbf{ngpts}}\left(l\right)$ is the total number of points on grid level $l$; $N={\mathbf{ngpts}}\left(l\right)×{\mathbf{npde}}$; $\Delta t$ is the current time step; ${u}_{t}$ is the time derivative of $u$ which is approximated by first-order finite differences; ${w}_{j}^{\mathrm{t}}$ is the time equivalent of the space weighting factor ${w}_{j}^{\mathrm{s}}$; and ${\alpha }_{ij}$ is given by
 $αij=τ ujmax100+uxi,yi,t$ (7)
where ${u}_{j}^{\mathrm{max}}$ is as before, and $\tau$ is the user-specified time tolerance.
An integration step is rejected and retried at all levels if
 $maxlμlt>1.0.$ (8)

## 10Example

For this routine two examples are presented, with a main program and two example problems given in Example 1 (EX1) and Example 2 (EX2).
Example 1 (EX1)
This example stems from combustion theory and is a model for a single, one-step reaction of a mixture of two chemicals (see Adjerid and Flaherty (1988)). The PDE for the temperature of the mixture $u$ is
 $∂u ∂t =d ∂2u ∂x2 + ∂2u ∂y2 +D1+α-u exp-δu$
for $0\le x,y\le 1$ and $t\ge 0$, with initial conditions $u\left(x,y,0\right)=1$ for $0\le x,y\le 1$, and boundary conditions
 $ux 0,y,t = 0, u 1,y,t = 1 for ​ 0≤y≤1,$
 $uy x,0,t = 0, u x,1,t = 1 for ​ 0≤x≤ 1.$
The heat release argument $\alpha =1$, the Damkohler number $D=R\mathrm{exp}\left(\delta \right)/\left(\alpha \delta \right)$, the activation energy $\delta =20$, the reaction rate $R=5$, and the diffusion argument $d=0.1$.
For small times the temperature gradually increases in a circular region about the origin, and at about $t=0.24$ ‘ignition’ occurs causing the temperature to suddenly jump from near unity to $1+\alpha$, and a reaction front forms and propagates outwards, becoming steeper. Thus during the solution, just one grid level is used up to the ignition point, then two levels, and then three as the reaction front steepens.
Example 2 (EX2)
This example is taken from a multispecies food web model, in which predator-prey relationships in a spatial domain are simulated (see Brown et al. (1994)). In this example there is just one species each of prey and predator, and the two PDEs for the concentrations ${c}_{1}$ and ${c}_{2}$ of the prey and the predator respectively are
 $∂c1 ∂t =c1b1+a11c1+a12c2+d1 ∂2c1 ∂x2 + ∂2c1 ∂y2 ,$
 $0=c2 b2+a21c1+a22c2+d2 ∂2 c2 ∂x2 + ∂2 c2 ∂y2 ,$
with
 $a11=a22=-1, a12=-0.5×10-6, and a21=104, and b1=1+αxy+βsin4πxsin4πy,$
where $\alpha =50$ and $\beta =300$, and ${b}_{2}=-{b}_{1}$.
The initial conditions are taken to be simple peaked functions which satisfy the boundary conditions and very nearly satisfy the PDEs:
 $c1=10+16x1-xy1-y2,$
 $c2=b2+a21c1,$
and the boundary conditions are of Neumann type, i.e., zero normal derivatives everywhere.
During the solution a number of peaks and troughs develop across the domain, and so the number of levels required increases with time. Since the solution varies rapidly in space across the whole of the domain, refinement at intermediate levels tends to occur at all points of the domain.

### 10.1Program Text

Program Text (d03rafe.f90)

### 10.2Program Data

Program Data (d03rafe.d)

### 10.3Program Results

Program Results (d03rafe.r)

© The Numerical Algorithms Group Ltd, Oxford, UK. 2017