NAG Library Routine Document
D03RBF
1 Purpose
D03RBF integrates a system of linear or nonlinear, timedependent partial differential equations (PDEs) in two space dimensions on a rectilinear 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 BiCGSTAB iterative linear solver with ILU preconditioning. Local uniform grid refinement is used to improve the accuracy of the solution. D03RBF originates from the VLUGR2 package (see
Blom and Verwer (1993) and
Blom et al. (1996)).
2 Specification
SUBROUTINE D03RBF ( 
NPDE, TS, TOUT, DT, TOLS, TOLT, INIDOM, PDEDEF, BNDARY, PDEIV, MONITR, OPTI, OPTR, RWK, LENRWK, IWK, LENIWK, LWK, LENLWK, ITRACE, IND, IFAIL) 
INTEGER 
NPDE, OPTI(4), LENRWK, IWK(LENIWK), LENIWK, LENLWK, ITRACE, IND, IFAIL 
REAL (KIND=nag_wp) 
TS, TOUT, DT(3), TOLS, TOLT, OPTR(3,NPDE), RWK(LENRWK) 
LOGICAL 
LWK(LENLWK) 
EXTERNAL 
INIDOM, PDEDEF, BNDARY, PDEIV, MONITR 

3 Description
D03RBF integrates the system of PDEs:
where
$\Omega $ is an arbitrary rectilinear domain, i.e., a domain bounded by perpendicular straight lines. If the domain is rectangular then it is recommended that
D03RAF is used.
The vector
$u$ is the set of solution values
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)$ for
$\left(x,y\right)\in \Omega $ must be specified at
$t={t}_{0}$ in
PDEIV.
Note that whilst complete generality is offered by the master equations
(1), D03RBF 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
The domain is covered by a uniform coarse base grid 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 usersupplied 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 8.
The system of PDEs and the boundary conditions are discretized in space on each grid using a standard secondorder finite difference scheme (centred on the internal domain and onesided at the boundaries), and the resulting system of ODEs is integrated in time using a secondorder, twostep, 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 usersupplied time tolerance
TOLT, and some further optional userspecified weighting factors held in
OPTR (see
Section 8 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) userspecified 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) userspecified 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 BiCGSTAB iterative method, in combination with ILU preconditioning. The maximum number of iterations can be specified by you in
OPTI.
In order to define the base grid you must first specify a virtual uniform rectangular grid which contains the entire base grid. The position of the virtual grid in physical $\left(x,y\right)$ space is given by the $\left(x,y\right)$ coordinates of its boundaries. The number of points ${n}_{x}$ and ${n}_{y}$ in the $x$ and $y$ directions must also be given, corresponding to the number of columns and rows respectively. This is sufficient to determine precisely the $\left(x,y\right)$ coordinates of all virtual grid points. Each virtual grid point is then referred to by integer coordinates $\left({v}_{x},{v}_{y}\right)$, where $\left(0,0\right)$ corresponds to the lowerleft corner and $\left({n}_{x}1,{n}_{y}1\right)$ corresponds to the upperright corner. ${v}_{x}$ and ${v}_{y}$ are also referred to as the virtual column and row indices respectively.
The base grid is then specified with respect to the virtual grid, with each base grid point coinciding with a virtual grid point. Each base grid point must be given an index, starting from $1$, and incrementing rowwise from the leftmost point of the lowest row. Also, each base grid row must be numbered consecutively from the lowest row in the grid, so that row $1$ contains grid point $1$.
As an example, consider the domain consisting of the two separate squares shown in
Figure 1. The lefthand diagram shows the virtual grid and its integer coordinates (i.e., its column and row indices), and the righthand diagram shows the base grid point indices and the base row indices (in brackets).
Hence the base grid point with index $6$ say is in base row $2$, virtual column 4, and virtual row $1$, i.e., virtual grid integer coordinates $\left(4,1\right)$; and the base grid point with index $19$ say is in base row $5$, virtual column 2, and virtual row $5$, i.e., virtual grid integer coordinates $\left(2,5\right)$.
The base grid must then be defined in
INIDOM by specifying the number of base grid rows, the number of base grid points, the number of boundaries, the number of boundary points, and the following integer arrays:
 LROW contains the base grid indices of the starting points of the base grid rows;
 IROW contains the virtual row numbers ${v}_{y}$ of the base grid rows;
 ICOL contains the virtual column numbers ${v}_{x}$ of the base grid points;
 LBND contains the grid indices of the boundary edges (without corners) and corner points;
 LLBND contains the starting elements of the boundaries and corners in LBND.
Finally,
ILBND contains the types of the boundaries and corners, as follows:
Boundaries:
 1 – lower boundary
 2 – left boundary
 3 – upper boundary
 4 – right boundary
External corners (
$90\xb0$):
 12 – lowerleft corner
 23 – upperleft corner
 34 – upperright corner
 41 – lowerright corner
Internal corners (
$270\xb0$):
 21 – lowerleft corner
 32 – upperleft corner
 43 – upperright corner
 14 – lowerright corner
Figure 2 shows the boundary types of a domain with a hole. Notice the logic behind the labelling of the corners: each one includes the types of the two adjacent boundary edges, in a clockwise fashion (outside the domain).
As an example, consider the domain shown in
Figure 3. The lefthand diagram shows the physical domain and the righthand diagram shows the base and virtual grids. The numbers outside the base grid are the indices of the left and rightmost base grid points, and the numbers inside the base grid are the boundary or corner numbers, indicating the order in which the boundaries are stored in
LBND.
For this example we have
NROWS = 11
NPTS = 105
NBNDS = 28
NBPTS = 72
LROW = (1,4,15,26,37,46,57,68,79,88,97)
IROW = (0,1,2,3,4,5,6,7,8,9,10)
ICOL = (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)
LBND = (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)
LLBND = (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)
ILBND = (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)
This particular domain is used in the example in
Section 9, and data statements are used to define the above arrays in that example program. For less complicated domains it is simpler to assign the values of the arrays in doloops. This also allows flexibility in the number of base grid points.
The routine
D03RYF can be called from
INIDOM to obtain a simple graphical representation of the base grid, and to verify the data that you have specified in
INIDOM.
Subgrids are stored internally using the same data structure, and solution information is communicated to you in
PDEIV,
PDEDEF and
BNDARY in arrays according to the grid index on the particular level, e.g.,
${\mathbf{X}}\left(i\right)$ and
${\mathbf{Y}}\left(i\right)$ contain the
$\left(x,y\right)$ coordinates of grid point
$i$, and
${\mathbf{U}}\left(i,j\right)$ contains the
$j$th solution component
${u}_{j}$ at grid point
$i$.
The grid data and the solutions at all grid levels are 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 provided. The usersupplied 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 parameter
TLAST. The routine
D03RZF should be called from
MONITR to obtain grid information at a particular level.
Further details of the underlying algorithm can be found in
Section 8 and in
Blom and Verwer (1993) and
Blom et al. (1996) and the references therein.
4 References
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 NMR9306 CWI, Amsterdam
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
5 Parameters
 1: NPDE – INTEGERInput
On entry: the number of PDEs in the system.
Constraint:
${\mathbf{NPDE}}\ge 1$.
 2: 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: TOUT – REAL (KIND=nag_wp)Input
On entry: the final value of $t$ to which the integration is to be carried out.
 4: DT($3$) – 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\times \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 ${\mathbf{DT}}\left(2\right)=10.0\times \mathit{machineprecision}$ 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$, $10.0\times \mathit{machineprecision}\times \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(\left{\mathbf{TS}}\right,\left{\mathbf{TOUT}}\right\right)\le {\mathbf{DT}}\left(1\right)\le {\mathbf{TOUT}}{\mathbf{TS}}$ 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: TOLS – REAL (KIND=nag_wp)Input
On entry: the space tolerance used in the grid refinement strategy (
$\sigma $ in equation
(4)). See
Section 8.2.
Constraint:
${\mathbf{TOLS}}>0.0$.
 6: TOLT – REAL (KIND=nag_wp)Input
On entry: the time tolerance used to determine the time step size (
$\tau $ in equation
(7)). See
Section 8.3.
Constraint:
${\mathbf{TOLT}}>0.0$.
 7: INIDOM – SUBROUTINE, supplied by the user.External Procedure
INIDOM must specify the base grid in terms of the data structure described in
Section 3.
INIDOM is not referenced if, on entry,
${\mathbf{IND}}=1$.
D03RYF can be called from
INIDOM to obtain a simple graphical representation of the base grid, and to verify the data that you have specified in
INIDOM. D03RBF also checks the validity of the data, but you are strongly advised to call
D03RYF to ensure that the base grid is exactly as required.
Note: the boundaries of the base grid should consist of as many points as are necessary to employ secondorder space discretization, i.e., a boundary enclosing the internal part of the domain must include at least $3$ grid points including the corners. If Neumann boundary conditions are to be applied the minimum is $4$.
The specification of
INIDOM is:
SUBROUTINE INIDOM ( 
MAXPTS, XMIN, XMAX, YMIN, YMAX, NX, NY, NPTS, NROWS, NBNDS, NBPTS, LROW, IROW, ICOL, LLBND, ILBND, LBND, IERR) 
INTEGER 
MAXPTS, NX, NY, NPTS, NROWS, NBNDS, NBPTS, LROW(*), IROW(*), ICOL(*), LLBND(*), ILBND(*), LBND(*), IERR 
REAL (KIND=nag_wp) 
XMIN, XMAX, YMIN, YMAX 

 1: MAXPTS – INTEGERInput
On entry: the maximum number of base grid points allowed by the available workspace.
 2: XMIN – REAL (KIND=nag_wp)Output
 3: XMAX – REAL (KIND=nag_wp)Output
On exit: the extents of the virtual grid 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.
 4: YMIN – REAL (KIND=nag_wp)Output
 5: YMAX – REAL (KIND=nag_wp)Output
On exit: the extents of the virtual grid in the $y$direction, i.e., the $y$ coordinates of the left and right boundaries respectively.
Constraint:
${\mathbf{YMIN}}<{\mathbf{YMAX}}$ and
YMAX must be sufficiently distinguishable from
YMIN for the precision of the machine being used.
 6: NX – INTEGEROutput
 7: NY – INTEGEROutput
On exit: the number of virtual grid points in the $x$ and $y$direction respectively (including the boundary points).
Constraint:
NX and
${\mathbf{NY}}\ge 4$.
 8: NPTS – INTEGEROutput
On exit: the total number of points in the base grid. If the required number of points is greater than
MAXPTS then
INIDOM must be exited immediately with
IERR set to
$1$ to avoid overwriting memory.
Constraint:
${\mathbf{NPTS}}\le {\mathbf{NX}}\times {\mathbf{NY}}$ and if ${\mathbf{IERR}}\ne 1$ on exit, ${\mathbf{NPTS}}\le {\mathbf{MAXPTS}}$.
 9: NROWS – INTEGEROutput
On exit: the total number of rows of the virtual grid that contain base grid points. This is the maximum base row index.
Constraint:
$4\le {\mathbf{NROWS}}\le {\mathbf{NY}}$.
 10: NBNDS – INTEGEROutput
On exit: the total number of physical boundaries and corners in the base grid.
Constraint:
${\mathbf{NBNDS}}\ge 8$.
 11: NBPTS – INTEGEROutput
On exit: the total number of boundary points in the base grid.
Constraint:
$12\le {\mathbf{NBPTS}}<{\mathbf{NPTS}}$.
 12: LROW($*$) – INTEGER arrayOutput
On exit: ${\mathbf{LROW}}\left(\mathit{i}\right)$, for $\mathit{i}=1,2,\dots ,{\mathbf{NROWS}}$, must contain the base grid index of the first grid point in base grid row $\mathit{i}$.
Constraints:
 $1\le {\mathbf{LROW}}\left(\mathit{i}\right)\le {\mathbf{NPTS}}$, for $\mathit{i}=1,2,\dots ,{\mathbf{NROWS}}$;
 ${\mathbf{LROW}}\left(\mathit{i}1\right)<{\mathbf{LROW}}\left(\mathit{i}\right)$, for $\mathit{i}=2,3,\dots ,{\mathbf{NROWS}}$.
 13: IROW($*$) – INTEGER arrayOutput
On exit: ${\mathbf{IROW}}\left(\mathit{i}\right)$, for $\mathit{i}=1,2,\dots ,{\mathbf{NROWS}}$, must contain the virtual row number ${v}_{y}$ that corresponds to base grid row $\mathit{i}$.
Constraints:
 $0\le {\mathbf{IROW}}\left(\mathit{i}\right)\le {\mathbf{NY}}$, for $\mathit{i}=1,2,\dots ,{\mathbf{NROWS}}$;
 ${\mathbf{IROW}}\left(\mathit{i}1\right)<{\mathbf{IROW}}\left(\mathit{i}\right)$, for $\mathit{i}=2,3,\dots ,{\mathbf{NROWS}}$.
 14: ICOL($*$) – INTEGER arrayOutput
On exit: ${\mathbf{ICOL}}\left(\mathit{i}\right)$, for $\mathit{i}=1,2,\dots ,{\mathbf{NPTS}}$, must contain the virtual column number ${v}_{x}$ that contains base grid point $\mathit{i}$.
Constraint:
$0\le {\mathbf{ICOL}}\left(\mathit{i}\right)\le {\mathbf{NX}}$, for $\mathit{i}=1,2,\dots ,{\mathbf{NPTS}}$.
 15: LLBND($*$) – INTEGER arrayOutput
On exit:
${\mathbf{LLBND}}\left(\mathit{i}\right)$, for
$\mathit{i}=1,2,\dots ,{\mathbf{NBNDS}}$, must contain the element of
LBND corresponding to the start of the
$\mathit{i}$th boundary or corner.
Note: the order of the boundaries and corners in
LLBND must be first all the boundaries and then all the corners. The end points of a boundary (i.e., the adjacent corner points) must
not be included in the list of points on that boundary. Also, if a corner is shared by two pairs of physical boundaries then it has two types and must therefore be treated as two corners.
Constraints:
 $1\le {\mathbf{LLBND}}\left(\mathit{i}\right)\le {\mathbf{NBPTS}}$, for $\mathit{i}=1,2,\dots ,{\mathbf{NBNDS}}$;
 ${\mathbf{LLBND}}\left(\mathit{i}1\right)<{\mathbf{LLBND}}\left(\mathit{i}\right)$, for $\mathit{i}=2,3,\dots ,{\mathbf{NBNDS}}$.
 16: ILBND($*$) – INTEGER arrayOutput
On exit:
${\mathbf{ILBND}}\left(\mathit{i}\right)$, for
$\mathit{i}=1,2,\dots ,{\mathbf{NBNDS}}$, must contain the type of the
$\mathit{i}$th boundary (or corner), as given in
Section 3.
Constraint:
${\mathbf{ILBND}}\left(\mathit{i}\right)$ must be equal to one of the following: $1$, $2$, $3$, $4$, $12$, $23$, $34$, $41$, $21$, $32$, $43$ or $14$, for $\mathit{i}=1,2,\dots ,{\mathbf{NBNDS}}$.
 17: LBND($*$) – INTEGER arrayOutput
On exit:
${\mathbf{LBND}}\left(\mathit{i}\right)$, for
$\mathit{i}=1,2,\dots ,{\mathbf{NBPTS}}$, must contain the grid index of the
$\mathit{i}$th boundary point. The order of the boundaries is as specified in
LLBND, but within this restriction the order of the points in
LBND is arbitrary.
Constraint:
$1\le {\mathbf{LBND}}\left(\mathit{i}\right)\le {\mathbf{NPTS}}$, for $\mathit{i}=1,2,\dots ,{\mathbf{NBPTS}}$.
 18: IERR – INTEGERInput/Output
On entry: will be initialized by D03RBF to some value prior to internal calls to
INIDOM.
On exit: if the required number of grid points is larger than
MAXPTS,
IERR must be set to
$1$ to force a termination of the integration and an immediate return to the calling program with
${\mathbf{IFAIL}}={\mathbf{3}}$. Otherwise,
IERR should remain unchanged.
INIDOM must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which D03RBF is called. Parameters denoted as
Input must
not be changed by this procedure.
 8: PDEDEF – SUBROUTINE, supplied by the user.External Procedure
PDEDEF must evaluate the functions
${F}_{\mathit{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:
SUBROUTINE PDEDEF ( 
NPTS, NPDE, T, X, Y, U, UT, UX, UY, UXX, UXY, UYY, RES) 
INTEGER 
NPTS, NPDE 
REAL (KIND=nag_wp) 
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), RES(NPTS,NPDE) 

 1: NPTS – INTEGERInput
On entry: the number of grid points in the current grid.
 2: NPDE – INTEGERInput
On entry: the number of PDEs in the system.
 3: T – REAL (KIND=nag_wp)Input
On entry: the current value of the independent variable $t$.
 4: X(NPTS) – 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: Y(NPTS) – 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: U(NPTS,NPDE) – 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: UT(NPTS,NPDE) – 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: UX(NPTS,NPDE) – 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: UY(NPTS,NPDE) – 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: UXX(NPTS,NPDE) – 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: UXY(NPTS,NPDE) – 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: UYY(NPTS,NPDE) – 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: RES(NPTS,NPDE) – 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 D03RBF is called. Parameters denoted as
Input must
not be changed by this procedure.
 9: 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:
SUBROUTINE BNDARY ( 
NPTS, NPDE, T, X, Y, U, UT, UX, UY, NBNDS, NBPTS, LLBND, ILBND, LBND, RES) 
INTEGER 
NPTS, NPDE, NBNDS, NBPTS, LLBND(NBNDS), ILBND(NBNDS), LBND(NBPTS) 
REAL (KIND=nag_wp) 
T, X(NPTS), Y(NPTS), U(NPTS,NPDE), UT(NPTS,NPDE), UX(NPTS,NPDE), UY(NPTS,NPDE), RES(NPTS,NPDE) 

 1: NPTS – INTEGERInput
On entry: the number of grid points in the current grid.
 2: NPDE – INTEGERInput
On entry: the number of PDEs in the system.
 3: T – REAL (KIND=nag_wp)Input
On entry: the current value of the independent variable $t$.
 4: X(NPTS) – 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: Y(NPTS) – 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: U(NPTS,NPDE) – 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: UT(NPTS,NPDE) – 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: UX(NPTS,NPDE) – 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: UY(NPTS,NPDE) – 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: NBNDS – INTEGERInput
On entry: the total number of physical boundaries and corners in the grid.
 11: NBPTS – INTEGERInput
On entry: the total number of boundary points in the grid.
 12: LLBND(NBNDS) – INTEGER arrayInput
On entry:
${\mathbf{LLBND}}\left(\mathit{i}\right)$, for
$\mathit{i}=1,2,\dots ,{\mathbf{NBNDS}}$, contains the element of
LBND corresponding to the start of the
$\mathit{i}$th boundary (or corner).
 13: ILBND(NBNDS) – INTEGER arrayInput
On entry:
${\mathbf{ILBND}}\left(\mathit{i}\right)$, for
$\mathit{i}=1,2,\dots ,{\mathbf{NBNDS}}$, contains the type of the
$\mathit{i}$th boundary, as given in
Section 3.
 14: LBND(NBPTS) – INTEGER arrayInput
On entry:
${\mathbf{LBND}}\left(\mathit{i}\right)$, contains the grid index of the
$\mathit{i}$th boundary point, where the order of the boundaries is as specified in
LLBND. 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),\mathit{j}\right)$, for
$\mathit{i}=1,2,\dots ,{\mathbf{NBPTS}}$ and
$\mathit{j}=1,2,\dots ,{\mathbf{NPDE}}$.
 15: RES(NPTS,NPDE) – REAL (KIND=nag_wp) arrayInput/Output
On entry: contains function values returned 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, i.e., points not included in
LBND, must
not be altered.
BNDARY must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which D03RBF is called. Parameters denoted as
Input must
not be changed by this procedure.
 10: PDEIV – SUBROUTINE, supplied by the user.External Procedure
PDEIV must specify the initial values of the PDE components
$u$ at all points in the base grid.
PDEIV is not referenced if, on entry,
${\mathbf{IND}}=1$.
The specification of
PDEIV is:
INTEGER 
NPTS, NPDE 
REAL (KIND=nag_wp) 
T, X(NPTS), Y(NPTS), U(NPTS,NPDE) 

 1: NPTS – INTEGERInput
On entry: the number of grid points in the base grid.
 2: NPDE – INTEGERInput
On entry: the number of PDEs in the system.
 3: T – REAL (KIND=nag_wp)Input
On entry: the (initial) value of the independent variable $t$.
 4: X(NPTS) – 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: Y(NPTS) – 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: U(NPTS,NPDE) – 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 D03RBF is called. Parameters denoted as
Input must
not be changed by this procedure.
 11: MONITR – SUBROUTINE, supplied by the user.External Procedure
MONITR is called by D03RBF 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 parameter
TLAST.
The input arguments contain information about the grid and solution at all grid levels used.
D03RZF should be called from
MONITR in order to extract the number of points and their
$\left(x,y\right)$ coordinates on a particular grid.
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:
SUBROUTINE MONITR ( 
NPDE, T, DT, DTNEW, TLAST, NLEV, XMIN, YMIN, DXB, DYB, LGRID, ISTRUC, LSOL, SOL, IERR) 
INTEGER 
NPDE, NLEV, LGRID(*), ISTRUC(*), LSOL(NLEV), IERR 
REAL (KIND=nag_wp) 
T, DT, DTNEW, XMIN, YMIN, DXB, DYB, SOL(*) 
LOGICAL 
TLAST 

 1: NPDE – INTEGERInput
On entry: the number of PDEs in the system.
 2: 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: 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: DTNEW – REAL (KIND=nag_wp)Input
On entry: the time step size that will be used for the next time step.
 5: 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: NLEV – INTEGERInput
On entry: the number of grid levels used at time
T.
 7: XMIN – REAL (KIND=nag_wp)Input
 8: YMIN – REAL (KIND=nag_wp)Input
On entry: the $\left(x,y\right)$ coordinates of the lowerleft corner of the virtual grid.
 9: DXB – REAL (KIND=nag_wp)Input
 10: DYB – REAL (KIND=nag_wp)Input
On entry: the sizes of the base grid spacing in the $x$ and $y$direction respectively.
 11: LGRID($*$) – INTEGER arrayInput
On entry: contains pointers to the start of the grid structures in
ISTRUC, and must be passed unchanged to
D03RZF in order to extract the grid information.
 12: ISTRUC($*$) – INTEGER arrayInput
On entry: contains the grid structures for each grid level and must be passed unchanged to
D03RZF in order to extract the grid information.
 13: LSOL(NLEV) – INTEGER arrayInput
On entry:
${\mathbf{LSOL}}\left(l\right)$ contains the pointer to the solution in
SOL at grid level
$l$ and time
T. (
${\mathbf{LSOL}}\left(l\right)$ actually contains the array index immediately preceding the start of the solution in
SOL.)
 14: SOL($*$) – REAL (KIND=nag_wp) arrayInput
On entry: contains the solution
$u$ at time
T for each grid level
$l$ in turn, positioned according to
LSOL.
More precisely
represents the
$\mathit{j}$th component of the solution at the
$\mathit{i}$th grid point in the
$\mathit{l}$th level, for
$\mathit{i}=1,2,\dots ,{n}_{\mathit{l}}$,
$\mathit{j}=1,2,\dots ,{\mathbf{NPDE}}$ and
$\mathit{l}=1,2,\dots ,{\mathbf{NLEV}}$, where
${n}_{\mathit{l}}$ is the number of grid points at level
$\mathit{l}$ (obtainable by a call to
D03RZF).
 15: IERR – INTEGERInput/Output
On entry: will be initialized by D03RBF to some value prior to internal calls to
IERR.
On exit: should be set to
$1$ to force a termination of the integration 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 D03RBF is called. Parameters denoted as
Input must
not be changed by this procedure.
 12: OPTI($4$) – 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$.
 13: OPTR($3$,NPDE) – REAL (KIND=nag_wp) arrayInput
On entry: may be used to specify the optional vectors
${u}^{\mathrm{max}}$,
${w}^{s}$ and
${w}^{t}$ in the space and time monitors (see
Section 8).
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}}^{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}}^{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}}$.
 14: RWK(LENRWK) – REAL (KIND=nag_wp) arrayCommunication Array
 15: LENRWK – INTEGERInput
On entry: the dimension of the array
RWK as declared in the (sub)program from which D03RBF is called.
The required value of
LENRWK cannot be determined exactly in advance, but a suggested value is
where
$l={\mathbf{OPTI}}\left(1\right)$ if
${\mathbf{OPTI}}\left(1\right)\ne 0$ and
$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).
Note: the size of
LENRWK cannot be checked upon initial entry to D03RBF since the number of grid points on the base grid is not known.
 16: IWK(LENIWK) – INTEGER arrayCommunication Array
On entry: if
${\mathbf{IND}}=0$,
IWK need not be set. Otherwise
IWK must remain unchanged from a previous call to D03RBF.
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
$l$ is a grid level taking the values
$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+l\right)$
 Contains the total number of residual evaluations performed (i.e., the number of times PDEDEF was called) at grid level $l$.
 ${\mathbf{IWK}}\left(2+m+l\right)$
 Contains the total number of Jacobian evaluations performed at grid level $l$.
 ${\mathbf{IWK}}\left(2+2\times m+l\right)$
 Contains the total number of Newton iterations performed at grid level $l$.
 ${\mathbf{IWK}}\left(2+3\times m+l\right)$
 Contains the total number of linear solver iterations performed at grid level $l$.
 ${\mathbf{IWK}}\left(2+4\times m+l\right)$
 Contains the maximum number of Newton iterations performed at any one time step at grid level $l$.
 ${\mathbf{IWK}}\left(2+5\times m+l\right)$
 Contains the maximum number of linear solver iterations performed at any one time step at grid level $l$.
Note: the total and maximum numbers are cumulative over all calls to D03RBF. If the specified maximum number of Newton or linear solver iterations is exceeded at any stage, then the maximums above are set to the specified maximum plus one.
 17: LENIWK – INTEGERInput
On entry: the dimension of the array
IWK as declared in the (sub)program from which D03RBF is called.
The required value of
LENIWK cannot be determined exactly in advance, but a suggested value is
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).
Note: the size of
LENIWK cannot be checked upon initial entry to D03RBF since the number of grid points on the base grid is not known.
 18: LWK(LENLWK) – LOGICAL arrayWorkspace
 19: LENLWK – INTEGERInput
On entry: the dimension of the array
LWK as declared in the (sub)program from which D03RBF is called.
The required value of
LENLWK cannot be determined exactly in advance, but a suggested value is
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).
Note: the size of
LENLWK cannot be checked upon initial entry to D03RBF since the number of grid points on the base grid is not known.
 20: ITRACE – INTEGERInput
On entry: the level of trace information required from D03RBF.
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$, then $1$ is assumed and similarly if ${\mathbf{ITRACE}}>3$, then $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.
 21: IND – INTEGERInput/Output
On entry: must be set to
$0$ or
$1$.
 ${\mathbf{IND}}=0$
 Starts the integration in time.
 ${\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 D03RBF: TOUT, ${\mathbf{DT}}\left(2\right)$, ${\mathbf{DT}}\left(3\right)$, TOLS, TOLT, OPTI, OPTR, ITRACE and IFAIL.
Constraint:
$0\le {\mathbf{IND}}\le 1$.
On exit: ${\mathbf{IND}}=1$.
Note: for users of the NAG Library for SMP & Multicore only, additional values of
IND are allowed. See
Section 2.2.3 in the Introduction to the NAG Library for SMP & Multicore for further details.
 22: IFAIL – INTEGERInput/Output

On entry:
IFAIL must be set to
$0$,
$1\text{ or}1$. If you are unfamiliar with this parameter you should refer to
Section 3.3 in the Essential Introduction 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 parameter, 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).
6 Error Indicators and Warnings
If on entry
${\mathbf{IFAIL}}={\mathbf{0}}$ or
${{\mathbf{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 $0.0<{\mathbf{DT}}\left(1\right)<10\times \mathit{machineprecision}\times \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(\left{\mathbf{TS}}\right,\left{\mathbf{TOUT}}\right\right)$, 
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)\text{ or}{\mathbf{DT}}\left(1\right)>{\mathbf{DT}}\left(3\right)$, 
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  ${\mathbf{IND}}\ne 0$ or $1$, 
or  ${\mathbf{IND}}=1$ on initial entry to D03RBF. 
 ${\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. At the initial time step this error may result because you set
IERR to
$1$ in
INIDOM or the internal check on the number of grid points following the call to
INIDOM. 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}}=6$
One or more of the output arguments of
INIDOM was incorrectly specified, i.e.,
 ${\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  ${\mathbf{NROWS}}<4$, 
or  ${\mathbf{NROWS}}>{\mathbf{NY}}$, 
or  ${\mathbf{NPTS}}>{\mathbf{NX}}\times {\mathbf{NY}}$, 
or  ${\mathbf{NBNDS}}<8$, 
or  ${\mathbf{NBPTS}}<12$, 
or  ${\mathbf{NBPTS}}\ge {\mathbf{NPTS}}$, 
or  ${\mathbf{LROW}}\left(i\right)<1$ or ${\mathbf{LROW}}\left(i\right)>{\mathbf{NPTS}}$, for some $i=1,2,\dots ,{\mathbf{NROWS}}$, 
or  ${\mathbf{LROW}}\left(i\right)\le {\mathbf{LROW}}\left(i1\right)$, for some $i=2,3,\dots ,{\mathbf{NROWS}}$, 
or  ${\mathbf{IROW}}\left(i\right)<0$ or ${\mathbf{IROW}}\left(i\right)>{\mathbf{NY}}$, for some $i=1,2,\dots ,{\mathbf{NROWS}}$, 
or  ${\mathbf{IROW}}\left(i\right)\le {\mathbf{IROW}}\left(i1\right)$, for some $i=2,3,\dots ,{\mathbf{NROWS}}$, 
or  ${\mathbf{ICOL}}\left(i\right)<0$ or ${\mathbf{ICOL}}\left(i\right)>{\mathbf{NX}}$, for some $i=1,2,\dots ,{\mathbf{NPTS}}$, 
or  ${\mathbf{LLBND}}\left(i\right)<1$ or ${\mathbf{LLBND}}\left(i\right)>{\mathbf{NBPTS}}$, for some $i=1,2,\dots ,{\mathbf{NBNDS}}$, 
or  ${\mathbf{LLBND}}\left(i\right)\le {\mathbf{LLBND}}\left(i1\right)$, for some $i=2,3,\dots ,{\mathbf{NBNDS}}$, 
or  ${\mathbf{ILBND}}\left(i\right)\ne 1$, $2$, $3$, $4$, $12$, $23$, $34$, $41$, $21$, $32$, $43$ or $14$, for some $i=1,2,\dots ,{\mathbf{NBNDS}}$, 
or  ${\mathbf{LBND}}\left(i\right)<1$ or ${\mathbf{LBND}}\left(i\right)>{\mathbf{NPTS}}$, for some $i=1,2,\dots ,{\mathbf{NBPTS}}$. 
7 Accuracy
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 parameters
TOLS and
TOLT described in
Section 8, and you should test the effects of varying these parameters. 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).
The local uniform grid refinement method is summarised as follows.
 Initialize the course base grid, an initial solution and an initial time step.
 Solve the system of PDEs on the current grid with the current time step.
 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$. 
 Update the coarser grid solution using the finer grid values.
 Estimate error in time integration. If time error is acceptable advance time level.
 Determine new step size then go to $2$ with coarse base as current grid.
For each grid point
$i$ a space monitor
${\mu}_{i}^{s}$ is determined by
where
$\Delta x$ and
$\Delta y$ are the grid widths in the
$x$ and
$y$ directions; and
${x}_{i}$,
${y}_{i}$ are the
$\left(x,y\right)$ coordinates at grid point
$i$. The parameter
${\gamma}_{j}$ is obtained from
where
$\sigma $ is the usersupplied 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}^{s}$ and
${u}_{j}^{\mathrm{max}}$ must also be supplied but may be set to the values
$1.0$ if little information about the solution is known.
A new level of refinement is created if
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 8.1.
The time integration is controlled using a time monitor calculated at each level
$l$ up to the maximum level used, given by
where
$\mathit{ngpts}\left(l\right)$ is the total number of points on grid level
$l$;
$N=\mathit{ngpts}\left(l\right)\times {\mathbf{NPDE}}$;
$\Delta t$ is the current time step;
${u}_{t}$ is the time derivative of
$u$ which is approximated by firstorder finite differences;
${w}_{j}^{t}$ is the time equivalent of the space weighting factor
${w}_{j}^{s}$; and
${\alpha}_{ij}$ is given by
where
${u}_{j}^{\mathrm{max}}$ is as before, and
$\tau $ is the userspecified time tolerance.
An integration step is rejected and retried at all levels if
9 Example
Note: a modified example is supplied with the NAG Library for SMP & Multicore and links have been supplied in the following subsections.
This example is taken from
Blom and Verwer (1993) and is the twodimensional Burgers' system
with
$\epsilon ={10}^{3}$ on the domain given in
Figure 3. Dirichlet boundary conditions are used on all boundaries using the exact solution
The solution contains a wave front at $y=x+0.25t$ which propagates in a direction perpendicular to the front with speed $\sqrt{2}/8$.
9.1 Program Text
Program Text (d03rbfe.f90)
Program Text (d03rbfe_smp.f90)
9.2 Program Data
Program Data (d03rbfe.d)
Program Data (d03rbfe_smp.d)
9.3 Program Results
Program Results (d03rbfe.r)
Program Results (d03rbfe_smp.r)