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

naginterfaces.library.pde.dim2_gen_order2_rectangle(ts, tout, dt, xmin, xmax, ymin, ymax, nx, ny, tols, tolt, pdedef, bndary, pdeiv, monitr, opti, optr, comm, itrace, ind, lenrwk=None, leniwk=None, lenlwk=None, data=None, io_manager=None, spiked_sorder='C')[source]

dim2_gen_order2_rectangle integrates a system of linear or nonlinear, time-dependent partial differential equations (PDEs) in two space dimensions on a rectangular domain. The method of lines is employed to reduce the PDEs to a system of ordinary differential equations (ODEs) which are solved using a backward differentiation formula (BDF) method. The resulting system of nonlinear equations is solved using a modified Newton method and a BI-CGSTAB iterative linear solver with ILU preconditioning. Local uniform grid refinement is used to improve the accuracy of the solution. dim2_gen_order2_rectangle 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 d03ra

https://www.nag.com/numeric/nl/nagdoc_27.3/flhtml/d03/d03raf.html

Parameters
tsfloat

The initial value of the independent variable .

toutfloat

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

dtfloat, array-like, shape

The initial, minimum and maximum time step sizes respectively.

Specifies the initial time step size to be used on the first entry, i.e., when . If then the default value is used. On subsequent entries (), the value of is not referenced.

Specifies the minimum time step size to be attempted by the integrator. If the default value is used.

Specifies the maximum time step size to be attempted by the integrator. If the default value is used.

xminfloat

The extents of the rectangular domain in the -direction, i.e., the coordinates of the left and right boundaries respectively.

xmaxfloat

The extents of the rectangular domain in the -direction, i.e., the coordinates of the left and right boundaries respectively.

yminfloat

The extents of the rectangular domain in the -direction, i.e., the coordinates of the lower and upper boundaries respectively.

ymaxfloat

The extents of the rectangular domain in the -direction, i.e., the coordinates of the lower and upper boundaries respectively.

nxint

The number of grid points in the -direction (including the boundary points).

nyint

The number of grid points in the -direction (including the boundary points).

tolsfloat

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

toltfloat

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

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

must evaluate the functions , for , 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 . is called for each subgrid in turn.

Parameters
tfloat

The current value of the independent variable .

xfloat, ndarray, shape

contains the coordinate of the th grid point, for .

yfloat, ndarray, shape

contains the coordinate of the th grid point, for .

ufloat, ndarray, shape

contains the value of the th PDE component at the th grid point, for , for .

utfloat, ndarray, shape

contains the value of for the th PDE component at the th grid point, for , for .

uxfloat, ndarray, shape

contains the value of for the th PDE component at the th grid point, for , for .

uyfloat, ndarray, shape

contains the value of for the th PDE component at the th grid point, for , for .

uxxfloat, ndarray, shape

contains the value of for the th PDE component at the th grid point, for , for .

uxyfloat, ndarray, shape

contains the value of for the th PDE component at the th grid point, for , for .

uyyfloat, ndarray, shape

contains the value of for the th PDE component at the th grid point, for , for .

dataarbitrary, optional, modifiable in place

User-communication data for callback functions.

Returns
resfloat, array-like, shape

must contain the value of , for , at the th grid point, for , 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, lbnd, res, data=None)

must evaluate the functions , for , 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 .

xfloat, ndarray, shape

contains the coordinate of the th grid point, for .

yfloat, ndarray, shape

contains the coordinate of the th grid point, for .

ufloat, ndarray, shape

contains the value of the th PDE component at the th grid point, for , for .

utfloat, ndarray, shape

contains the value of for the th PDE component at the th grid point, for , for .

uxfloat, ndarray, shape

contains the value of for the th PDE component at the th grid point, for , for .

uyfloat, ndarray, shape

contains the value of for the th PDE component at the th grid point, for , for .

lbndint, ndarray, shape

contains the grid index for the th boundary point, for . Hence the th boundary point has coordinates and , and the corresponding solution values are , etc.

resfloat, ndarray, shape

contains the value of , for , at the th grid point, for , as returned by . The residuals at the boundary points will be overwritten and so need not have been set by .

dataarbitrary, optional, modifiable in place

User-communication data for callback functions.

Returns
resfloat, array-like, shape

must contain the value of , for , at the th boundary point, for .

Note: elements of corresponding to interior points must not be altered.

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

must specify the initial values of the PDE components at all points in the grid. is not referenced if, on entry, .

Parameters
npdeint

The number of PDEs in the system.

tfloat

The (initial) value of the independent variable .

xfloat, ndarray, shape

contains the coordinate of the th grid point, for .

yfloat, ndarray, shape

contains the coordinate of the th grid point, for .

dataarbitrary, optional, modifiable in place

User-communication data for callback functions.

Returns
ufloat, array-like, shape

must contain the value of the th PDE component at the th grid point, for , for .

monitrcallable ierr = monitr(npde, t, dt, dtnew, tlast, ngpts, xpts, ypts, lsol, sol, ierr, data=None)

is called by dim2_gen_order2_rectangle 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 .

The input arguments contain information about the grid and solution at all grid levels used.

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 , i.e., the time at the end of the integration step just completed.

dtfloat

The current time step size , i.e., the time step size used for the integration step just completed.

dtnewfloat

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

tlastbool

Indicates if intermediate or final time step. for an intermediate step, for the last call to before returning to your program.

ngptsint, ndarray, shape

contains the number of grid points at level , for .

xptsfloat, ndarray, shape

Contains the coordinates of the grid points in each level in turn, i.e., , for , for .

So for level , , where , for , for .

yptsfloat, ndarray, shape

Contains the coordinates of the grid points in each level in turn, i.e., , for , for .

So for level , , where , for , for .

lsolint, ndarray, shape

contains the pointer to the solution in at grid level and time . ( actually contains the array index immediately preceding the start of the solution in .)

solfloat, ndarray, shape

Contains the solution at time for each grid level in turn, positioned according to , i.e., for level , , for , for , for .

ierrint

Will be set to .

dataarbitrary, optional, modifiable in place

User-communication data for callback functions.

Returns
ierrint

Should be set to to force a tidy termination and an immediate return to the calling program with = 4. should remain unchanged otherwise.

optiint, array-like, shape

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

All the default options are employed.

The default value of , for , can be obtained by setting .

Specifies the maximum number of grid levels allowed (including the base grid). . The default value is .

Specifies the maximum number of Jacobian evaluations allowed during each nonlinear equations solution. . The default value is .

Specifies the maximum number of Newton iterations in each nonlinear equations solution. . The default value is .

Specifies the maximum number of iterations in each linear equations solution. . The default value is .

optrfloat, array-like, shape

May be used to specify the optional vectors , and in the space and time monitors (see Further Comments).

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

, for , specifies , the approximate maximum absolute value of the th component of , as used in (4) and (7). , for .

, for , specifies , the weighting factors used in the space monitor (see (4)) to indicate the relative importance of the th component of on the space monitor. , for .

, for , specifies , the weighting factors used in the time monitor (see (6)) to indicate the relative importance of the th component of on the time monitor. , for .

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_rectangle. may take the value , , , or .

No output is generated.

Only warning messages are printed.

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

If , is assumed and similarly if , is assumed.

The advisory messages are given in greater detail as increases.

Setting allows you to monitor the progress of the integration without possibly excessive information.

indint

Must be set to or , alternatively or .

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

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: , , , , , , and . is assumed to be serial.

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

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: , , , , , , and . is assumed to have been parallelized by you, as described in Parallelism and Performance. In all other respects, this is equivalent to .

lenrwkNone or int, optional

Note: if this argument is None then a default value will be used, determined as follows: .

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

where if and otherwise, and 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 = 3 and an estimated required size is printed.

leniwkNone or int, optional

Note: if this argument is None then a default value will be used, determined as follows: .

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

where is the expected maximum number of grid points at any one level and if and otherwise.

If during the execution the supplied value is found to be too small then the function returns with = 3 and an estimated required size is printed.

lenlwkNone or int, optional

Note: if this argument is None then a default value will be used, determined as follows: .

The required value of cannot be determined exactly in advanced, but a suggested value is

where 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 = 3 and an estimated required size is printed.

dataarbitrary, optional

User-communication data for callback functions.

io_managerFileObjManager, optional

Manager for I/O in this routine.

spiked_sorderstr, optional

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

spiked_sorder =

row-major storage will be used;

spiked_sorder =

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 which has been reached. Normally .

dtfloat, ndarray, shape

contains the time step size for the next time step. and are unchanged or set to their default values if zero on entry.

indint

, if on input was or , or , if on input was or .

Raises
NagValueError
(errno )

On entry, is not equal to or : .

(errno )

On entry, too small for initial grid level: , minimum value . Note that subsequent levels will require more. Consult document.

(errno )

On entry, too small for initial grid level: , minimum value . Note that subsequent levels will require more. Consult document.

(errno )

On entry, too small for initial grid level: , minimum value . Note that subsequent levels will require more. Consult document.

(errno )

On entry, too close to : and .

(errno )

On entry, too close to : and .

(errno )

On entry, and .

Constraint: .

(errno )

On entry, and .

Constraint: .

(errno )

On entry, .

Constraint: .

(errno )

On entry, .

Constraint: .

(errno )

On entry, .

Constraint: .

(errno )

On entry, and . Note that was reset to default if zero on entry.

Constraint: if , .

(errno )

On entry, and . Note that was reset to default if zero on entry.

Constraint: if , .

(errno )

On entry, and too large: and .

(errno )

On entry, and too small: and .

(errno )

On entry, .

Constraint: if , .

(errno )

On entry, and .

Constraint: .

(errno )

On entry, .

Constraint: .

(errno )

On entry, .

Constraint: .

(errno )

On entry, , and .

Constraint: .

(errno )

On entry, and .

Constraint: .

(errno )

On entry, , and .

Constraint: if , .

(errno )

On entry, .

Constraint: .

(errno )

On entry, .

Constraint: .

(errno )

On entry, .

Constraint: .

(errno )

On entry, too small: .

(errno )

On entry, and .

Constraint: .

(errno )

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

(errno )

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

Warns
NagAlgorithmicWarning
(errno )

set to in . Integration completed as far as : .

(errno )

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_rectangle integrates the system of PDEs:

for and in the rectangular domain , , and time interval , where the vector is the set of solution values

and denotes partial differentiation with respect to , and similarly for etc.

The functions must be supplied by you in . Similarly the initial values of the functions must be specified at in .

Note that whilst complete generality is offered by the master equations (1), dim2_gen_order2_rectangle is not appropriate for all PDEs. In particular, hyperbolic systems should not be solved using this function. Also, at least one component of must appear in the system of PDEs.

The boundary conditions must be supplied by you in in the form

for all when or and for all when or and

The domain is covered by a uniform coarse base grid of size 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 . 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 and . 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 , and some further optional user-specified weighting factors held in (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 . The step size is adjusted so that the remaining integration interval is an integer number times . In this way a solution is obtained at .

A modified Newton method is used to solve the nonlinear equations arising from the time integration. You may specify (in ) 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 . 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 .

The solution at all grid levels is stored in the workspace arrays, along with other information needed for a restart (i.e., a continuation call). It is not intended that you extract the solution from these arrays, indeed the necessary information regarding these arrays is not included. The user-supplied monitor should be used to obtain the solution at particular levels and times. is called at the end of every time step, with the last step being identified via the input argument .

Within , , and the data structure is as follows. Each point on a particular grid is given an index (ranging from to the total number of points on the grid) and all coordinate or solution information is stored in arrays according to this index, e.g., and contain the - and coordinate of point , and contains the th solution component at point .

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

Adjerid, S and Flaherty, J E, 1988, A local refinement finite element method for two-dimensional parabolic systems, SIAM J. Sci. Statist. Comput. (9), 792–811

Blom, J G, Trompert, R A and Verwer, J G, 1996, Algorithm 758. VLUGR2: A vectorizable adaptive grid solver for PDEs in 2D, Trans. Math. Software (22), 302–328

Blom, J G and Verwer, J G, 1993, VLUGR2: A vectorized local uniform grid refinement code for PDEs in 2D, Report NM-R9306, CWI, Amsterdam

Brown, P N, Hindmarsh, A C and Petzold, L R, 1994, Using Krylov methods in the solution of large scale differential-algebraic systems, SIAM J. Sci. Statist. Comput. (15), 1467–1488

Trompert, R A, 1993, Local uniform grid refinement and systems of coupled partial differential equations, Appl. Numer. Maths (12), 331–355

Trompert, R A and Verwer, J G, 1993, Analysis of the implicit Euler local uniform grid refinement method, SIAM J. Sci. Comput. (14), 259–278