naginterfaces.library.pde.dim2_​gen_​order2_​rectilinear

naginterfaces.library.pde.dim2_gen_order2_rectilinear(ts, tout, dt, tols, tolt, inidom, pdedef, bndary, pdeiv, monitr, opti, optr, comm, itrace, ind, lenrwk=None, leniwk=None, lenlwk=None, data=None, io_manager=None, spiked_sorder='C')

dim2_gen_order2_rectilinear integrates a system of linear or nonlinear, time-dependent 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 BI-CGSTAB iterative linear solver with ILU preconditioning. Local uniform grid refinement is used to improve the accuracy of the solution. dim2_gen_order2_rectilinear originates from the VLUGR2 package (see Blom and Verwer (1993) and Blom et al. (1996)).

For full information please refer to the NAG Library document for d03rb

https://support.nag.com/numeric/nl/nagdoc_30.2/flhtml/d03/d03rbf.html

Parameters

tsfloat

The initial value of the independent variable $t$.

toutfloat

The final value of $t$ to which the integration is to be carried out.

dtfloat, array-like, shape $\left(3\right)$

The initial, minimum and maximum time step sizes respectively.

$dt\left[0\right]$

Specifies the initial time step size to be used on the first entry, i.e., when $ind=0$. If $dt\left[0\right]=0.0$ then the default value $dt\left[0\right]=0.01\times (tout-ts)$ is used. On subsequent entries ($ind=1$), the value of $dt\left[0\right]$ is not referenced.

$dt\left[1\right]$

Specifies the minimum time step size to be attempted by the integrator. If $dt\left[1\right]=0.0$ the default value $dt\left[1\right]=10.0\times \text{machine precision}$ is used.

$dt\left[2\right]$

Specifies the maximum time step size to be attempted by the integrator. If $dt\left[2\right]=0.0$ the default value $dt\left[2\right]=tout-ts$ is used.

tolsfloat

The space tolerance used in the grid refinement strategy ($\sigma$ in equation (4)). See Further Comments.

toltfloat

The time tolerance used to determine the time step size ($\tau$ in equation (7)). See Time Integration.

inidomcallable (xmin, xmax, ymin, ymax, nx, ny, npts, nrows, nbnds, nbpts, lrow, irow, icol, llbnd, ilbnd, lbnd, ierr) = inidom(maxpts, ierr, data=None)

$inidom$ must specify the base grid in terms of the data structure described in Notes. $inidom$ is not referenced if, on entry, $ind=1$.

Note: the boundaries of the base grid should consist of as many points as are necessary to employ second-order 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$.

Parameters

maxptsint

The maximum number of base grid points allowed by the available workspace.

ierrint

Will be initialized by dim2_gen_order2_rectilinear to some value prior to internal calls to $inidom$.

dataarbitrary, optional, modifiable in place

User-communication data for callback functions.

Returns

xminfloat

The extents of the virtual grid in the $x$-direction, i.e., the $x$ coordinates of the left and right boundaries respectively.

xmaxfloat

The extents of the virtual grid in the $x$-direction, i.e., the $x$ coordinates of the left and right boundaries respectively.

yminfloat

The extents of the virtual grid in the $y$-direction, i.e., the $y$ coordinates of the left and right boundaries respectively.

ymaxfloat

The extents of the virtual grid in the $y$-direction, i.e., the $y$ coordinates of the left and right boundaries respectively.

nxint

The number of virtual grid points in the $x$- and $y$-direction respectively (including the boundary points).

nyint

The number of virtual grid points in the $x$- and $y$-direction respectively (including the boundary points).

nptsint

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.

nrowsint

The total number of rows of the virtual grid that contain base grid points. This is the maximum base row index.

nbndsint

The total number of physical boundaries and corners in the base grid.

nbptsint

The total number of boundary points in the base grid.

lrowint, array-like, shape $\left(maxpts\right)$

$lrow[\text{i}-1]$, for $\text{i}=1,2,\dots ,nrows$, must contain the base grid index of the first grid point in base grid row $\text{i}$.

irowint, array-like, shape $\left(maxpts\right)$

$irow[\text{i}-1]$, for $\text{i}=1,2,\dots ,nrows$, must contain the virtual row number $v_y$ that corresponds to base grid row $\text{i}$.

icolint, array-like, shape $\left(maxpts\right)$

$icol[\text{i}-1]$, for $\text{i}=1,2,\dots ,npts$, must contain the virtual column number $v_x$ that contains base grid point $\text{i}$.

llbndint, array-like, shape $\left(maxpts\right)$

$llbnd[\text{i}-1]$, for $\text{i}=1,2,\dots ,nbnds$, must contain the element of $lbnd$ corresponding to the start of the $\text{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.

ilbndint, array-like, shape $\left(maxpts\right)$

$ilbnd[\text{i}-1]$, for $\text{i}=1,2,\dots ,nbnds$, must contain the type of the $\text{i}$th boundary (or corner), as given in Notes.

lbndint, array-like, shape $\left(maxpts\right)$

$lbnd[\text{i}-1]$, for $\text{i}=1,2,\dots ,nbpts$, must contain the grid index of the $\text{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.

ierrint

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 $errno$ = 3. Otherwise, $ierr$ should remain unchanged.

pdedefcallable res = pdedef(t, x, y, u, ut, ux, uy, uxx, uxy, uyy, data=None)

$pdedef$ must evaluate the functions $F_{\text{j}}$, for $\text{j}=1,2,\dots ,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.

Parameters

tfloat

The current value of the independent variable $t$.

xfloat, ndarray, shape $\left(\text{npts}\right)$

$x[\text{i}-1]$ contains the $x$ coordinate of the $\text{i}$th grid point, for $\text{i}=1,2,\dots ,npts$.

yfloat, ndarray, shape $\left(\text{npts}\right)$

$y[\text{i}-1]$ contains the $y$ coordinate of the $\text{i}$th grid point, for $\text{i}=1,2,\dots ,npts$.

ufloat, ndarray, shape $(\text{npts},\text{npde})$

$u[\text{i}-1,\text{j}-1]$ contains the value of the $\text{j}$th PDE component at the $\text{i}$th grid point, for $\text{j}=1,2,\dots ,npde$, for $\text{i}=1,2,\dots ,npts$.

utfloat, ndarray, shape $(\text{npts},\text{npde})$

$ut[\text{i}-1,\text{j}-1]$ contains the value of $\frac{\partial u}{\partial t}$ for the $\text{j}$th PDE component at the $\text{i}$th grid point, for $\text{j}=1,2,\dots ,npde$, for $\text{i}=1,2,\dots ,npts$.

uxfloat, ndarray, shape $(\text{npts},\text{npde})$

$ux[\text{i}-1,\text{j}-1]$ contains the value of $\frac{\partial u}{\partial x}$ for the $\text{j}$th PDE component at the $\text{i}$th grid point, for $\text{j}=1,2,\dots ,npde$, for $\text{i}=1,2,\dots ,npts$.

uyfloat, ndarray, shape $(\text{npts},\text{npde})$

$uy[\text{i}-1,\text{j}-1]$ contains the value of $\frac{\partial u}{\partial y}$ for the $\text{j}$th PDE component at the $\text{i}$th grid point, for $\text{j}=1,2,\dots ,npde$, for $\text{i}=1,2,\dots ,npts$.

uxxfloat, ndarray, shape $(\text{npts},\text{npde})$

$uxx[\text{i}-1,\text{j}-1]$ contains the value of $\frac{{\partial }^{2}u}{\partial x^2}$ for the $\text{j}$th PDE component at the $\text{i}$th grid point, for $\text{j}=1,2,\dots ,npde$, for $\text{i}=1,2,\dots ,npts$.

uxyfloat, ndarray, shape $(\text{npts},\text{npde})$

$uxy[\text{i}-1,\text{j}-1]$ contains the value of $\frac{{\partial }^{2}u}{\partial x\partial y}$ for the $\text{j}$th PDE component at the $\text{i}$th grid point, for $\text{j}=1,2,\dots ,npde$, for $\text{i}=1,2,\dots ,npts$.

uyyfloat, ndarray, shape $(\text{npts},\text{npde})$

$uyy[\text{i}-1,\text{j}-1]$ contains the value of $\frac{{\partial }^{2}u}{\partial y^2}$ for the $\text{j}$th PDE component at the $\text{i}$th grid point, for $\text{j}=1,2,\dots ,npde$, for $\text{i}=1,2,\dots ,npts$.

dataarbitrary, optional, modifiable in place

User-communication data for callback functions.

Returns

resfloat, array-like, shape $(\text{npts},\text{npde})$

$res[\text{i}-1,\text{j}-1]$ must contain the value of $F_{\text{j}}$, for $\text{j}=1,2,\dots ,npde$, at the $\text{i}$th grid point, for $\text{i}=1,2,\dots ,npts$, although the residuals at boundary points will be ignored (and overwritten later on) and so they need not be specified here.

bndarycallable res = bndary(t, x, y, u, ut, ux, uy, llbnd, ilbnd, lbnd, res, data=None)

$bndary$ must evaluate the functions $G_{\text{j}}$, for $\text{j}=1,2,\dots ,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 function.

Parameters

tfloat

The current value of the independent variable $t$.

xfloat, ndarray, shape $\left(\text{npts}\right)$

$x[\text{i}-1]$ contains the $x$ coordinate of the $\text{i}$th grid point, for $\text{i}=1,2,\dots ,npts$.

yfloat, ndarray, shape $\left(\text{npts}\right)$

$y[\text{i}-1]$ contains the $y$ coordinate of the $\text{i}$th grid point, for $\text{i}=1,2,\dots ,npts$.

ufloat, ndarray, shape $(\text{npts},\text{npde})$

$u[\text{i}-1,\text{j}-1]$ contains the value of the $\text{j}$th PDE component at the $\text{i}$th grid point, for $\text{j}=1,2,\dots ,npde$, for $\text{i}=1,2,\dots ,npts$.

utfloat, ndarray, shape $(\text{npts},\text{npde})$

$ut[\text{i}-1,\text{j}-1]$ contains the value of $\frac{\partial u}{\partial t}$ for the $\text{j}$th PDE component at the $\text{i}$th grid point, for $\text{j}=1,2,\dots ,npde$, for $\text{i}=1,2,\dots ,npts$.

uxfloat, ndarray, shape $(\text{npts},\text{npde})$

$ux[\text{i}-1,\text{j}-1]$ contains the value of $\frac{\partial u}{\partial x}$ for the $\text{j}$th PDE component at the $\text{i}$th grid point, for $\text{j}=1,2,\dots ,npde$, for $\text{i}=1,2,\dots ,npts$.

uyfloat, ndarray, shape $(\text{npts},\text{npde})$

$uy[\text{i}-1,\text{j}-1]$ contains the value of $\frac{\partial u}{\partial y}$ for the $\text{j}$th PDE component at the $\text{i}$th grid point, for $\text{j}=1,2,\dots ,npde$, for $\text{i}=1,2,\dots ,npts$.

llbndint, ndarray, shape $\left(\text{nbnds}\right)$

$llbnd[\text{i}-1]$, for $\text{i}=1,2,\dots ,nbnds$, contains the element of $lbnd$ corresponding to the start of the $\text{i}$th boundary (or corner).

ilbndint, ndarray, shape $\left(\text{nbnds}\right)$

$ilbnd[\text{i}-1]$, for $\text{i}=1,2,\dots ,nbnds$, contains the type of the $\text{i}$th boundary, as given in Notes.

lbndint, ndarray, shape $\left(\text{nbpts}\right)$

$lbnd[\text{i}-1]$, contains the grid index of the $\text{i}$th boundary point, where the order of the boundaries is as specified in $llbnd$. Hence the $\text{i}$th boundary point has coordinates $x\left[lbnd\right[\text{i}-1]-1]$ and $y\left[lbnd\right[\text{i}-1]-1]$, and the corresponding solution values are $u\left[lbnd\right[\text{i}-1]-1,\text{j}-1]$, for $\text{j}=1,2,\dots ,npde$, for $\text{i}=1,2,\dots ,nbpts$.

resfloat, ndarray, shape $(\text{npts},\text{npde})$

Contains function values returned by $pdedef$.

dataarbitrary, optional, modifiable in place

User-communication data for callback functions.

Returns

resfloat, array-like, shape $(\text{npts},\text{npde})$

$res\left[lbnd\right[\text{i}-1]-1,\text{j}-1]$ must contain the value of $G_{\text{j}}$, for $\text{j}=1,2,\dots ,npde$, at the $\text{i}$th boundary point, for $\text{i}=1,2,\dots ,nbpts$.

Note: elements of $res$ corresponding to interior points, i.e., points not included in $lbnd$, must not be altered.

pdeivcallable u = pdeiv(npde, t, x, y, data=None)

$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, $ind=1$.

Parameters

npdeint

The number of PDEs in the system.

tfloat

The (initial) value of the independent variable $t$.

xfloat, ndarray, shape $\left(\text{npts}\right)$

$x[\text{i}-1]$ contains the $x$ coordinate of the $\text{i}$th grid point, for $\text{i}=1,2,\dots ,npts$.

yfloat, ndarray, shape $\left(\text{npts}\right)$

$y[\text{i}-1]$ contains the $y$ coordinate of the $\text{i}$th grid point, for $\text{i}=1,2,\dots ,npts$.

dataarbitrary, optional, modifiable in place

User-communication data for callback functions.

Returns

ufloat, array-like, shape $(\text{npts},npde)$

$u[\text{i}-1,\text{j}-1]$ must contain the value of the $\text{j}$th PDE component at the $\text{i}$th grid point, for $\text{j}=1,2,\dots ,npde$, for $\text{i}=1,2,\dots ,npts$.

monitrcallable ierr = monitr(npde, t, dt, dtnew, tlast, xmin, ymin, dxb, dyb, lgrid, istruc, lsol, sol, ierr, data=None)

$monitr$ is called by dim2_gen_order2_rectilinear 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. dim2_gen_order2_rectilinear_extractgrid() should be called from $monitr$ in order to extract the number of points and their $(x,y)$ 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.

Parameters

npdeint

The number of PDEs in the system.

tfloat

The current value of the independent variable $t$, i.e., the time at the end of the integration step just completed.

dtfloat

The current time step size $\Delta t$, i.e., the time step size used for the integration step just completed.

dtnewfloat

The time step size that will be used for the next time step.

tlastbool

Indicates if intermediate or final time step. $tlast=False$ for an intermediate step, $tlast=True$ for the last call to $monitr$ before returning to your program.

xminfloat

The $(x,y)$ coordinates of the lower-left corner of the virtual grid.

yminfloat

The $(x,y)$ coordinates of the lower-left corner of the virtual grid.

dxbfloat

The sizes of the base grid spacing in the $x$- and $y$-direction respectively.

dybfloat

The sizes of the base grid spacing in the $x$- and $y$-direction respectively.

lgridint, ndarray, shape $(\text{nlev}+1)$

Contains pointers to the start of the grid structures in $istruc$, and must be passed unchanged to dim2_gen_order2_rectilinear_extractgrid() in order to extract the grid information.

istrucint, ndarray, shape $\left(\text{nistruc}\right)$

Contains the grid structures for each grid level and must be passed unchanged to dim2_gen_order2_rectilinear_extractgrid() in order to extract the grid information.

lsolint, ndarray, shape $\left(\text{nlev}\right)$

$lsol[l-1]$ contains the pointer to the solution in $sol$ at grid level $l$ and time $t$. ($lsol[l-1]$ actually contains the array index immediately preceding the start of the solution in $sol$.)

solfloat, ndarray, shape $\left(\text{nsol}\right)$

Contains the solution $u$ at time $t$ for each grid level $l$ in turn, positioned according to $lsol$. More precisely

$$u[\text{i}-1,\text{j}-1]=sol\left[lsol\right[\text{l}-1]+(\text{j}-1)\times n_{\text{l}}+\text{i}-1]$$

represents the $\text{j}$th component of the solution at the $\text{i}$th grid point in the $\text{l}$th level, for $\text{l}=1,2,\dots ,\text{nlev}$, for $\text{j}=1,2,\dots ,npde$, for $\text{i}=1,2,\dots ,n_{\text{l}}$, where $n_{\text{l}}$ is the number of grid points at level $\text{l}$ (obtainable by a call to dim2_gen_order2_rectilinear_extractgrid()).

ierrint

Will be initialized by dim2_gen_order2_rectilinear to some value prior to internal calls to $ierr$.

dataarbitrary, optional, modifiable in place

User-communication data for callback functions.

Returns

ierrint

Should be set to $1$ to force a termination of the integration and an immediate return to the calling program with $errno$ = 4. $ierr$ should remain unchanged otherwise.

optiint, array-like, shape $\left(4\right)$

May be set to control various options available in the integrator.

$opti\left[0\right]=0$

All the default options are employed.

$opti\left[0\right]>0$

The default value of $opti[\text{i}-1]$, for $\text{i}=2,3,\dots ,4$, can be obtained by setting $opti[\text{i}-1]=0$.

$opti\left[0\right]$

Specifies the maximum number of grid levels allowed (including the base grid). $opti\left[0\right]\ge 0$. The default value is $opti\left[0\right]=3$.

$opti\left[1\right]$

Specifies the maximum number of Jacobian evaluations allowed during each nonlinear equations solution. $opti\left[1\right]\ge 0$. The default value is $opti\left[1\right]=2$.

$opti\left[2\right]$

Specifies the maximum number of Newton iterations in each nonlinear equations solution. $opti\left[2\right]\ge 0$. The default value is $opti\left[2\right]=10$.

$opti\left[3\right]$

Specifies the maximum number of iterations in each linear equations solution. $opti\left[3\right]\ge 0$. The default value is $opti\left[3\right]=100$.

optrfloat, array-like, shape $(3,\text{npde})$

May be used to specify the optional vectors $u^{max}$, $w^s$ and $w^t$ in the space and time monitors (see Further Comments).

If an optional vector is not required then all its components should be set to $1.0$.

$optr[0,\text{j}-1]$, for $\text{j}=1,2,\dots ,\text{npde}$, specifies $u_{\text{j}}^{max}$, the approximate maximum absolute value of the $\text{j}$th component of $u$, as used in (4) and (7). $optr[0,\text{j}-1]>0.0$, for $\text{j}=1,2,\dots ,\text{npde}$.

$optr[1,\text{j}-1]$, for $\text{j}=1,2,\dots ,\text{npde}$, specifies $w_{\text{j}}^s$, the weighting factors used in the space monitor (see (4)) to indicate the relative importance of the $\text{j}$th component of $u$ on the space monitor. $optr[1,\text{j}-1]\ge 0.0$, for $\text{j}=1,2,\dots ,\text{npde}$.

$optr[2,\text{j}-1]$, for $\text{j}=1,2,\dots ,\text{npde}$, specifies $w_{\text{j}}^t$, the weighting factors used in the time monitor (see (6)) to indicate the relative importance of the $\text{j}$th component of $u$ on the time monitor. $optr[2,\text{j}-1]\ge 0.0$, for $\text{j}=1,2,\dots ,\text{npde}$.

commdict, communication object, modified in place

Communication structure.

On initial entry: need not be set.

itraceint

The level of trace information required from dim2_gen_order2_rectilinear. $itrace$ may take the value $-1$, $0$, $1$, $2$ or $3$.

$itrace=-1$

No output is generated.

$itrace=0$

Only warning messages are printed.

$itrace>0$

Output from the underlying solver is printed. This output contains details of the time integration, the nonlinear iteration and the linear solver.

If $itrace<-1$, $-1$ is assumed and similarly if $itrace>3$, $3$ is assumed.

The advisory messages are given in greater detail as $itrace$ increases.

Setting $itrace=1$ allows you to monitor the progress of the integration without possibly excessive information.

indint

Must be set to $0$ or $1$, alternatively $10$ or $11$.

$ind=0$

Starts the integration in time. $pdedef$ is assumed to be serial.

$ind=1$

Continues the integration after an earlier exit from the function. In this case, only the following parameters may be reset between calls to dim2_gen_order2_rectilinear: $\text{tout}$, $\text{dt}$, $\text{tols}$, $\text{tolt}$, $\text{opti}$, $\text{optr}$, $\text{itrace}$ and $\text{errno}$. $pdedef$ is assumed to be serial.

$ind=10$

Starts the integration in time. $pdedef$ is assumed to have been parallelized by you, as described in Parallelism and Performance. In all other respects, this is equivalent to $ind=0$.

$ind=11$

Continues the integration after an earlier exit from the function. In this case, only the following parameters may be reset between calls to dim2_gen_order2_rectangle(): $\text{tout}$, $\text{dt}$, $\text{tols}$, $\text{tolt}$, $\text{opti}$, $\text{optr}$, $\text{itrace}$ and $\text{errno}$. $pdedef$ is assumed to have been parallelized by you, as described in Parallelism and Performance. In all other respects, this is equivalent to $ind=1$.

lenrwkNone or int, optional

Note: if this argument is None then a default value will be used, determined as follows: $\text{maxpts}\times npde\times (5\times \text{maxlev}+18\times npde+9)+2\times \text{maxpts}$.

The required value of $lenrwk$ cannot be determined exactly in advance, but a suggested value is

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

where $l=opti\left[0\right]$ if $opti\left[0\right]\ne 0$ and $l=3$ otherwise, and $\text{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 function returns with $errno$ = 3 and an estimated required size is printed on the current error message unit (see FileObjManager).

Note: the size of $lenrwk$ cannot be checked upon initial entry to dim2_gen_order2_rectilinear since the number of grid points on the base grid is not known.

leniwkNone or int, optional

Note: if this argument is None then a default value will be used, determined as follows: $\text{maxpts}\times (14+5\times \text{maxlev})+7+2$.

The required value of $leniwk$ cannot be determined exactly in advance, but a suggested value is

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

where $\text{maxpts}$ is the expected maximum number of grid points at any one level and $m=opti\left[0\right]$ if $opti\left[0\right]>0$ and $m=3$ otherwise.

If during the execution the supplied value is found to be too small then the function returns with $errno$ = 3 and an estimated required size is printed on the current error message unit (see FileObjManager).

Note: the size of $leniwk$ cannot be checked upon initial entry to dim2_gen_order2_rectilinear since the number of grid points on the base grid is not known.

lenlwkNone or int, optional

Note: if this argument is None then a default value will be used, determined as follows: $\text{maxpts}+1$.

The required value of $lenlwk$ cannot be determined exactly in advance, but a suggested value is

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

where $\text{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 function returns with $errno$ = 3 and an estimated required size is printed.

Note: the size of $lenlwk$ cannot be checked upon initial entry to dim2_gen_order2_rectilinear since the number of grid points on the base grid is not known.

dataarbitrary, optional

User-communication data for callback functions.

io_managerFileObjManager, optional

Manager for I/O in this routine.

spiked_sorderstr, optional

If $optr$ is spiked (i.e., has unit extent in all but one dimension, or has size $1$), $spiked\_sorder$ selects the storage order to associate with it in the NAG Engine:

spiked_sorder = $\text{'C'}$

row-major storage will be used;

spiked_sorder = $\text{'F'}$

column-major storage will be used.

Two-dimensional arrays returned from callback functions in this routine must then use the same storage order.

Returns

tsfloat

The value of $t$ which has been reached. Normally $ts=tout$.

dtfloat, ndarray, shape $\left(3\right)$

$dt\left[0\right]$ contains the time step size for the next time step. $dt\left[1\right]$ and $dt\left[2\right]$ are unchanged or set to their default values if zero on entry.

indint

$ind=1$, if $ind$ on input was $0$ or $1$, or $ind=11$, if $ind$ on input was $10$ or $11$.

Raises

NagValueError

(errno $1$)

On entry, $ind$ is not equal to $0$ or $1$: $ind=⟨value⟩$.

(errno $1$)

On entry, $\text{npde}=⟨value⟩$.

Constraint: $\text{npde}\ge 1$.

(errno $1$)

On entry, $dt\left[0\right]=⟨value⟩$ and $dt\left[2\right]=⟨value⟩$. Note that $dt\left[2\right]$ was reset to default if zero on entry.

Constraint: if $ind=0$, $dt\left[0\right]\le dt\left[2\right]$.

(errno $1$)

On entry, $dt\left[0\right]=⟨value⟩$ and $dt\left[1\right]=⟨value⟩$. Note that $dt\left[1\right]$ was reset to default if zero on entry.

Constraint: if $ind=0$, $dt\left[0\right]\ge dt\left[1\right]$.

(errno $1$)

On entry, $ind=0$ and $dt\left[0\right]$ too large: $dt\left[0\right]=⟨value⟩$ and $\text{maximum value}=⟨value⟩$.

(errno $1$)

On entry, $ind=0$ and $dt\left[0\right]$ too small: $dt\left[0\right]=⟨value⟩$ and $\text{minimum value}=⟨value⟩$.

(errno $1$)

On entry, $dt\left[0\right]=⟨value⟩$.

Constraint: if $ind=0$, $dt\left[0\right]\ge 0.0$.

(errno $1$)

On entry, $dt\left[1\right]=⟨value⟩$ and $dt\left[2\right]=⟨value⟩$.

Constraint: $dt\left[1\right]\le dt\left[2\right]$.

(errno $1$)

On entry, $dt\left[2\right]=⟨value⟩$.

Constraint: $dt\left[2\right]\ge 0.0$.

(errno $1$)

On entry, $dt\left[1\right]=⟨value⟩$.

Constraint: $dt\left[1\right]\ge 0.0$.

(errno $1$)

On entry, $\text{i}=⟨value⟩$, $\text{j}=⟨value⟩$ and $optr[\text{i}-1,\text{j}-1]=⟨value⟩$.

Constraint: $optr[\text{i}-1,\text{j}-1]\ge 0.0$.

(errno $1$)

On entry, $\text{j}=⟨value⟩$ and $optr[0,\text{j}-1]=⟨value⟩$.

Constraint: $optr[0,\text{j}-1]>0.0$.

(errno $1$)

On entry, $opti\left[0\right]=⟨value⟩$, $\text{i}=⟨value⟩$ and $opti[\text{i}-1]=⟨value⟩$.

Constraint: if $opti\left[0\right]>0$, $opti[\text{i}-1]\ge 0$.

(errno $1$)

On entry, $opti\left[0\right]=⟨value⟩$.

Constraint: $opti\left[0\right]\ge 0$.

(errno $1$)

On entry, $tolt=⟨value⟩$.

Constraint: $tolt>0.0$.

(errno $1$)

On entry, $tols=⟨value⟩$.

Constraint: $tols>0.0$.

(errno $1$)

On entry, $tout-ts$ too small: $tout-ts=⟨value⟩$.

(errno $1$)

On entry, $tout=⟨value⟩$ and $ts=⟨value⟩$.

Constraint: $tout>ts$.

(errno $2$)

Attempted time-step smaller than specified minimum. Check problem formulation in $pdedef$, $bndary$ and $pdeiv$. Try increasing $itrace$ for more information.

(errno $3$)

One or more of the workspace arrays are too small. Try increasing $itrace$ for more information.

(errno $6$)

Invalid output argument from $inidom$.

Warns

NagAlgorithmicWarning

(errno $4$)

$IERR$ set to $1$ in $monitr$. Integration completed as far as $ts$: $ts=⟨value⟩$.

(errno $5$)

Integration completed, but maximum number of levels too small for required accuracy.

Notes

No equivalent traditional C interface for this routine exists in the NAG Library.

dim2_gen_order2_rectilinear integrates the system of PDEs:

$$F_j(t,x,y,u,u_t,u_x,u_y,u_{xx},u_{xy},u_{yy})=0\text{,}j=1,2,\dots ,\text{npde}\text{,}(x,y)\in \Omega \text{,}{t}_0\le t\le t_{out}\text{,}$$

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 dim2_gen_order2_rectangle() is used.

The vector $u$ is the set of solution values

$$u(x,y,t)={[u_1(x,y,t),\dots ,u_{\text{npde}}(x,y,t)]}^T\text{,}$$

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(x,y,t)$ for $(x,y)\in \Omega$ must be specified at $t=t_0$ in $pdeiv$.

Note that whilst complete generality is offered by the master equations (1), dim2_gen_order2_rectilinear is not appropriate for all PDEs. In particular, hyperbolic systems should not be solved using this function. 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

$$G_j(t,x,y,u,u_t,u_x,u_y)=0\text{,}j=1,2,\dots ,\text{npde}\text{,}(x,y)\in \partial \Omega \text{,}{t}_0\le t\le t_{out}\text{.}$$

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 user-supplied space tolerance $tols$. A number of options, 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 Further Comments.

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 Further Comments 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_{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$.

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 $(x,y)$ space is given by the $(x,y)$ 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 $(x,y)$ coordinates of all virtual grid points. Each virtual grid point is then referred to by integer coordinates $(v_x,v_y)$, where $(0,0)$ corresponds to the lower-left corner and $(n_x-1,n_y-1)$ corresponds to the upper-right 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 row-wise 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 [label omitted]. The left-hand diagram shows the virtual grid and its integer coordinates (i.e., its column and row indices), and the right-hand diagram shows the base grid point indices and the base row indices (in brackets).

[figure omitted]

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 $(4,1)$; 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 $(2,5)$.

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°$):

12 – lower-left corner

23 – upper-left corner

34 – upper-right corner

41 – lower-right corner

Internal corners ($270°$):

21 – lower-left corner

32 – upper-left corner

43 – upper-right corner

14 – lower-right corner

Figure [label omitted] 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).

[figure omitted]

As an example, consider the domain shown in Figure [label omitted]. The left-hand diagram shows the physical domain and the right-hand 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$.

[figure omitted]

For this example we would

set $nrows$ to $11$

set $npts$ to $105$

set $nbnds$ to $28$

set $nbpts$ to $72$

set $lrow$ to the vector $(1,4,15,26,37,46,57,68,79,88,97)$

set $irow$ to the vector $(0,1,2,3,4,5,6,7,8,9,10)$

set $icol$ to the vector $(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)$

set $lbnd$ to the vector $(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)$

set $llbnd$ to the vector $(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)$

set $ilbnd$ to the vector $(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)$

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., $x[i-1]$ and $y[i-1]$ contain the $(x,y)$ coordinates of grid point $i$, and $u[i-1,j-1]$ 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 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$. The function dim2_gen_order2_rectilinear_extractgrid() should be called from $monitr$ to obtain grid information at a particular level.

Further details of the underlying algorithm can be found in Further Comments and in Blom and Verwer (1993) and Blom et al. (1996) and the references therein.

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 NM-R9306, 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