Linear Partial Differential Equations

The LinearPDE class is used to define a general linear, steady, second order PDE for an unknown function $ u$ on a given $ \Omega$ defined through a Domain object. In the following $ \Gamma$ denotes the boundary of the domain $ \Omega$. $ n$ denotes the outer normal field on $ \Gamma$.

For a single PDE with a solution with a single component the linear PDE is defined in the following form:

$\displaystyle -(A\hackscore{jl} u\hackscore{,l})\hackscore{,j}-(B\hackscore{j} u)\hackscore{,j}+C\hackscore{l} u\hackscore{,l}+D u =-X\hackscore{j,j}+Y \; .$ (65)

$ u_{,j}$ denotes the derivative of $ u$ with respect to the $ j$-th spatial direction. Einstein's summation convention, ie. summation over indexes appearing twice in a term of a sum is performed, is used. The coefficients $ A$, $ B$, $ C$, $ D$, $ X$ and $ Y$ have to be specified through Data objects in the general FunctionSpace on the PDE or objects that can be converted into such Data objects. $ A$ is a rank-2 Data object, $ B$, $ C$ and $ X$ are rank-1 Data object and $ D$ and $ Y$ are scalar. The following natural boundary conditions are considered on $ \Gamma$:

$\displaystyle n\hackscore{j}(A\hackscore{jl} u\hackscore{,l}+B\hackscore{j} u)+d u=n\hackscore{j}X\hackscore{j} + y \;.$ (66)

Notice that the coefficients $ A$, $ B$ and $ X$ are defined in the PDE. The coefficients $ d$ and $ y$ are each a scalar Data object in the boundary FunctionSpace. Constraints for the solution prescribing the value of the solution at certain locations in the domain. They have the form

$\displaystyle u=r$    where $\displaystyle q>0$ (67)

$ r$ and $ q$ are each scalar Data object where $ q$ is the characteristic function defining where the constraint is applied. The constraints defined by Equation (4.3) override any other condition set by Equation (4.1) or Equation (4.2).

For a system of PDEs and a solution with several components the PDE has the form

$\displaystyle -(A\hackscore{ijkl} u\hackscore{k,l})\hackscore{,j}-(B\hackscore{...
...ore{k,l}+D\hackscore{ik} u\hackscore{k} =-X\hackscore{ij,j}+Y\hackscore{i} \; .$ (68)

$ A$ is a rank-4 Data object, $ B$ and $ C$ are each a rank-3 Data object, $ D$ and $ X$ are each a rank-2 Data object and $ Y$ is a rank-1 Data object. The natural boundary conditions take the form:

$\displaystyle n\hackscore{j}(A\hackscore{ijkl} u\hackscore{k,l}+B\hackscore{ijk...
...d\hackscore{ik} u\hackscore{k}=n\hackscore{j}X\hackscore{ij}+y\hackscore{i} \;.$ (69)

The coefficient $ d$ is a rank-2 Data object and $ y$ is a rank-1 Data object both in the boundary FunctionSpace. Constraints take the form

$\displaystyle u\hackscore{i}=r\hackscore{i}$    where $\displaystyle q\hackscore{i}>0$ (70)

$ r$ and $ q$ are each rank-1 Data object. Notice that not necessarily all components must have a constraint at all locations.

LinearPDE also supports solution discontinuities over contact region $ \Gamma^{contact}$ in the domain $ \Omega$. To specify the conditions across the discontinuity we are using the generalised flux $ J$[*] which is in the case of a systems of PDEs and several components of the solution defined as

$\displaystyle J\hackscore{ij}=A\hackscore{ijkl}u\hackscore{k,l}+B\hackscore{ijk}u\hackscore{k}-X\hackscore{ij}$ (71)

For the case of single solution component and single PDE $ J$ is defined

$\displaystyle J\hackscore{j}=A\hackscore{jl}u\hackscore{,l}+B\hackscore{j}u\hackscore{k}-X\hackscore{j}$ (72)

In the context of discontinuities $ n$ denotes the normal on the discontinuity pointing from side 0 towards side 1. For a system of PDEs the contact condition takes the form

$\displaystyle n\hackscore{j} J^{0}\hackscore{ij}=n\hackscore{j} J^{1}\hackscore{ij}=y^{contact}\hackscore{i} - d^{contact}\hackscore{ik} [u]\hackscore{k} \; .$ (73)

where $ J^{0}$ and $ J^{1}$ are the fluxes on side 0 and side $ 1$ of the discontinuity $ \Gamma^{contact}$, respectively. $ [u]$, which is the difference of the solution at side 1 and at side 0, denotes the jump of $ u$ across $ \Gamma^{contact}$. The coefficient $ d^{contact}$ is a rank-2 Data object and $ y^{contact}$ is a rank-1 Data object both in the contact FunctionSpace on side 0 or contact FunctionSpace on side 1. In case of a single PDE and a single component solution the contact condition takes the form

$\displaystyle n\hackscore{j} J^{0}\hackscore{j}=n\hackscore{j} J^{1}\hackscore{j}=y^{contact} - d^{contact}[u]$ (74)

In this case the the coefficient $ d^{contact}$ and $ y^{contact}$ are each scalar Data object both in the contact FunctionSpace on side 0 or contact FunctionSpace on side 1.

The PDE is symmetrical if

$\displaystyle A\hackscore{jl}=A\hackscore{lj}$    and $\displaystyle B\hackscore{j}=C\hackscore{j}$ (75)

The system of PDEs is symmetrical if
$\displaystyle A\hackscore{ijkl}$ $\displaystyle =$ $\displaystyle A\hackscore{klij}$ (76)
$\displaystyle B\hackscore{ijk}$ $\displaystyle =$ $\displaystyle C\hackscore{kij}$ (77)
$\displaystyle D\hackscore{ik}$ $\displaystyle =$ $\displaystyle D\hackscore{ki}$ (78)
$\displaystyle d\hackscore{ik}$ $\displaystyle =$ $\displaystyle d\hackscore{ki}$ (79)
$\displaystyle d^{contact}\hackscore{ik}$ $\displaystyle =$ $\displaystyle d^{contact}\hackscore{ki}$ (80)

Note that in contrast with the scalar case Equation (4.11) now the coefficients $ D$, $ d$ abd $ d^{contact}$ have to be inspected.

The following example illustrates the typical usage of the LinearPDE class:
\begin{python}
from esys.escript import *
from esys.escript.linearPDEs import Li...
...setValue(A=kappa*kronecker(mydomain),D=1,Y=1)
u=mypde.getSolution()
\end{python}
We refer to chapter 1 for more details.

An instance of the SolverOptions class is attached to the LinearPDE class object. It is used to set options of the solver used to solve the PDE. In the following code the getSolverOptions is used to access the SolverOptions attached to mypde:
\begin{python}
from esys.escript import *
from esys.escript.linearPDEs import Li...
...-8)
mypde.getSolverOptions().setIterMax(1000)
u=mypde.getSolution()
\end{python}
In this code the preconditoned conjugate gradient method SolverOptions.PCG with preconditioner SolverOptions.AMG . The relative tolerance is set tto $ 10^{-8}$ and the maximum number of iteration steps to $ 1000$.

Moreover, after a completed solution call the attached SolverOptions object gives access to diagnostic informations:
\begin{python}
u=mypde.getSolution()
print 'Number of iteration steps =', mypde....
...norm of returned solution =', mypde.getDiagnostics('residual_norm')
\end{python}
Typically a negative value for a diagnostic value indicates that the value is undefined.



Subsections
esys@esscc.uq.edu.au