NAG Library Routine Document

d03raf  (dim2_gen_order2_rectangle)

1

Purpose

2

Specification

3

Description

4

References

5

Arguments

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\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{machine precision}$ 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{machine precision}\times \max \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:     $\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)

C Header Interface

#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, nbpts, 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)

C Header Interface

#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)

C Header Interface

#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, dtnew, tlast, nlev, ngpts, 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

C Header Interface

#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)\times \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^{\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}}^{\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

$$\mathbf{lenrwk}=\mathit{maxpts}\times \mathbf{npde}\times \left(5\times \mathit{l}+18\times \mathbf{npde}+9\right)+2\times \mathit{maxpts}\text{,}$$

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}\times \mathbf{ny}\times \mathbf{npde}\times \left(14+18\times \mathbf{npde}\right)+2\times \mathbf{nx}\times \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\times m+\mathit{l}\right)$

Contains the total number of Newton iterations performed at grid level $\mathit{l}$.

$\mathbf{iwk}\left(2+3\times m+\mathit{l}\right)$

Contains the total number of linear solver iterations performed at grid level $\mathit{l}$.

$\mathbf{iwk}\left(2+4\times 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\times 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

$$\mathbf{leniwk}=\mathit{maxpts}\times \left(14+5\times m\right)+7\times m+2\text{,}$$

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\times \mathbf{nx}\times \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

$$\mathbf{lenlwk}=\mathit{maxpts}+1\text{,}$$

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}\times \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).

6

Error Indicators and Warnings

$\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{machine precision}\times \max \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)$ 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$

An unexpected error has been triggered by this routine. Please contact NAG.

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.

7

Accuracy

8

Parallelism and Performance

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))

!$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

9

Further Comments

9.1

Algorithm Outline

9.2

Refinement Strategy

9.3

Time Integration

10

Example

10.1

Program Text

Program Text (d03rafe.f90)

10.2

Program Data

Program Data (d03rafe.d)

10.3

Program Results

Program Results (d03rafe.r)