naginterfaces.library.pde.dim2_​ellip_​mgrid

naginterfaces.library.pde.dim2_ellip_mgrid(ngx, ngy, a, rhs, ub, maxit, acc, iout, io_manager=None)

dim2_ellip_mgrid solves seven-diagonal systems of linear equations which arise from the discretization of an elliptic partial differential equation on a rectangular region. This function uses a multigrid technique.

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

https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/d03/d03edf.html

Parameters

ngxint

The number of interior grid points in the $x$-direction, $n_x$. $ngx-1$ should preferably be divisible by as high a power of $2$ as possible.

ngyint

The number of interior grid points in the $y$-direction, $n_y$. $ngy-1$ should preferably be divisible by as high a power of $2$ as possible.

afloat, array-like, shape $((4\times (ngx+1)\times (ngy+1))/3,7)$

$a[\text{i}+(\text{j}-1)\times ngx-1,\text{k}-1]$ must be set to $a_{\text{i}\text{j}}^{\text{k}}$, for $\text{k}=1,2,\dots ,7$, for $\text{j}=1,2,\dots ,ngy$, for $\text{i}=1,2,\dots ,ngx$.

rhsfloat, array-like, shape $\left((4\times (ngx+1)\times (ngy+1))/3\right)$

$rhs[\text{i}+(\text{j}-1)\times ngx-1]$ must be set to $f_{\text{i}\text{j}}$, for $\text{j}=1,2,\dots ,ngy$, for $\text{i}=1,2,\dots ,ngx$.

ubfloat, array-like, shape $(ngx\times ngy)$

$ub[i+(j-1)\times ngx-1]$ must be set to the initial estimate for the solution $u_{ij}$.

maxitint

The maximum permitted number of multigrid iterations. If $maxit=0$, no multigrid iterations are performed, but the coarse-grid approximations and incomplete Crout decompositions are computed, and may be output if $iout$ is set accordingly.

accfloat

The required tolerance for convergence of the residual $2$-norm:

$${\parallel r\parallel }_2=\sqrt{\sum _{k=1}^{ngx\times ngy}{\left(r_k\right)}^2}$$

where $r=f-Au$ and $u$ is the computed solution. Note that the norm is not scaled by the number of equations. The function will stop after fewer than $maxit$ iterations if the residual $2$-norm is less than the specified tolerance. (If $maxit>0$, at least one iteration is always performed.)

If on entry $acc=0.0$, the machine precision is used as a default value for the tolerance; if $acc>0.0$, but $acc$ is less than the machine precision, the function will stop when the residual $2$-norm is less than the machine precision and $\text{errno}$ will be set to $4$.

ioutint

Controls the output of printed information to the file object associated with the advisory I/O unit (see FileObjManager):

$iout=0$

No output.

$iout=1$

The solution $u_{\text{i}\text{j}}$, for $\text{j}=1,2,\dots ,ngy$, for $\text{i}=1,2,\dots ,ngx$.

$iout=2$

The residual $2$-norm after each iteration, with the reduction factor over the previous iteration.

$iout=3$

As for $iout=1$ and $iout=2$.

$iout=4$

As for $iout=3$, plus the final residual (as returned in $ub$).

$iout=5$

As for $iout=4$, plus the initial elements of $a$ and $rhs$.

$iout=6$

As for $iout=5$, plus the Galerkin coarse grid approximations.

$iout=7$

As for $iout=6$, plus the incomplete Crout decompositions.

$iout=8$

As for $iout=7$, plus the residual after each iteration.

The elements $a[p-1,k-1]$, the Galerkin coarse grid approximations and the incomplete Crout decompositions are output in the format:

Y-index $=j$

X-index $=ia[p-1,0]a[p-1,1]a[p-1,2]a[p-1,3]a[p-1,4]a[p-1,5]a[p-1,6]$

where $p=\text{i}+(\text{j}-1)\times ngx$, for $\text{j}=1,2,\dots ,ngy$, for $\text{i}=1,2,\dots ,ngx$.

The vectors $u[p-1]$, $ub[p-1]$, $rhs[p-1]$ are output in matrix form with $ngy$ rows and $ngx$ columns.

Where $ngx>10$, the $ngx$ values for a given $j$ value are produced in rows of $10$.

Values of $iout>4$ may yield considerable amounts of output.

io_managerFileObjManager, optional

Manager for I/O in this routine.

Returns

afloat, ndarray, shape $((4\times (ngx+1)\times (ngy+1))/3,7)$

Is overwritten.

rhsfloat, ndarray, shape $\left((4\times (ngx+1)\times (ngy+1))/3\right)$

The first $ngx\times ngy$ elements are unchanged and the rest of the array is used as workspace.

ubfloat, ndarray, shape $(ngx\times ngy)$

The corresponding component of the residual $r=f-au$.

usfloat, ndarray, shape $\left((4\times (ngx+1)\times (ngy+1))/3\right)$

The residual $2$-norm, stored in element $us\left[0\right]$.

ufloat, ndarray, shape $\left((4\times (ngx+1)\times (ngy+1))/3\right)$

The computed solution $u_{\text{i}\text{j}}$ is returned in $u[\text{i}+(\text{j}-1)\times ngx-1]$, for $\text{j}=1,2,\dots ,ngy$, for $\text{i}=1,2,\dots ,ngx$.

numitint

The number of iterations performed.

Raises

NagValueError

(errno $1$)

On entry, $maxit=⟨value⟩$.

Constraint: $maxit\ge 0$.

(errno $1$)

On entry, $acc=⟨value⟩$.

Constraint: $acc\ge 0.0$.

(errno $1$)

On entry, $iout=⟨value⟩$.

Constraint: $0\le iout\le 8$.

(errno $1$)

On entry, $ngy=⟨value⟩$.

Constraint: $ngy\ge 3$.

(errno $1$)

On entry, $ngx=⟨value⟩$.

Constraint: $ngx\ge 3$.

(errno $3$)

After $maxit$ iterations the residual norm is not less than the tolerance $maxit=⟨value⟩$, residual norm $=⟨value⟩$, tolerance $=⟨value⟩$. The residual norm increased at one or more iterations after the first.

Warns

NagAlgorithmicWarning

(errno $2$)

After $maxit$ iterations the residual norm is not less than the tolerance $maxit=⟨value⟩$, residual norm $=⟨value⟩$, tolerance $=⟨value⟩$. The residual norm has decreased at each iteration after the first.

(errno $4$)

On entry, $acc$ is less than machine precision. The function terminated because the residual norm is less than machine precision. residual norm $=⟨value⟩$, machine precision $=⟨value⟩$ and $acc=⟨value⟩$.

Notes

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

dim2_ellip_mgrid solves, by multigrid iteration, the seven-point scheme

$$\begin{array}{c}\begin{array}{cccccc}& A_{i,j}^6{u}_{i-1,j+1}& +& A_{i,j}^7{u}_{i,j+1}& & \\ +& A_{i,j}^3{u}_{i-1,j}& +& A_{i,j}^4{u}_{ij}& +& A_{i,j}^5{u}_{i+1,j}\\ & & +& A_{i,j}^1{u}_{i,j-1}& +& A_{i,j}^2{u}_{i+1,j-1}=f_{ij}\text{,}i=1,2,\dots ,n_x\text{and}j=1,2,\dots ,n_y\text{,}\end{array}\end{array}$$

which arises from the discretization of an elliptic partial differential equation of the form

$$\alpha (x,y)U_{xx}+\beta (x,y)U_{xy}+\gamma (x,y)U_{yy}+\delta (x,y)U_x+ϵ(x,y)U_y+\varphi (x,y)U=\psi (x,y)$$

and its boundary conditions, defined on a rectangular region. This we write in matrix form as

$$Au=f\text{.}$$

The algorithm is described in separate reports by Wesseling (1982a), Wesseling (1982b) and McCarthy (1983).

Systems of linear equations, matching the seven-point stencil defined above, are solved by a multigrid iteration. An initial estimate of the solution must be provided by you. A zero guess may be supplied if no better approximation is available.

A ‘smoother’ based on incomplete Crout decomposition is used to eliminate the high frequency components of the error. A restriction operator is then used to map the system on to a sequence of coarser grids. The errors are then smoothed and prolongated (mapped onto successively finer grids). When the finest cycle is reached, the approximation to the solution is corrected. The cycle is repeated for $maxit$ iterations or until the required accuracy, $acc$, is reached.

dim2_ellip_mgrid will automatically determine the number $l$ of possible coarse grids, ‘levels’ of the multigrid scheme, for a particular problem. In other words, dim2_ellip_mgrid determines the maximum integer $l$ so that $n_x$ and $n_y$ can be expressed in the form

$$n_x=m{2}^{l-1}+1\text{,}{n}_y=n{2}^{l-1}+1\text{, with}m\ge 2\text{and}n\ge 2\text{.}$$

It should be noted that the rate of convergence improves significantly with the number of levels used (see McCarthy (1983)), so that $n_x$ and $n_y$ should be carefully chosen so that $n_x-1$ and $n_y-1$ have factors of the form $2^l$, with $l$ as large as possible. For good convergence the integer $l$ should be at least $2$.

dim2_ellip_mgrid has been found to be robust in application, but being an iterative method the problem of divergence can arise. For a strictly diagonally dominant matrix $A$

$$\left|A_{ij}^4\right|>\sum _{k\ne 4}\left|A_{ij}^k\right|\text{,}i=1,2,\dots ,n_x\text{and}j=1,2,\dots ,n_y$$

no such problem is foreseen. The diagonal dominance of $A$ is not a necessary condition, but should this condition be strongly violated then divergence may occur. The quickest test is to try the function.

References

McCarthy, G J, 1983, Investigation into the multigrid code MGD1, Report AERE-R 10889, Harwell

Wesseling, P, 1982, MGD1 – a robust and efficient multigrid method, Multigrid Methods. Lecture Notes in Mathematics (960), 614–630, Springer–Verlag

Wesseling, P, 1982, Theoretical aspects of a multigrid method, SIAM J. Sci. Statist. Comput. (3), 387–407