doc/user/linearPDE.tex

% $Id$
%
%           Copyright © 2006 by ACcESS MNRF
%               \url{http://www.access.edu.au
%         Primary Business: Queensland, Australia.
%   Licensed under the Open Software License version 3.0
%      http://www.opensource.org/licenses/osl-3.0.php
%


\chapter{The module \linearPDEs}

\declaremodule{extension}{linearPDEs} \modulesynopsis{Linear partial pifferential equation handler}
The module \linearPDEs provides an interface to define and solve linear partial 
differential equations within \escript. \linearPDEs does not provide any 
solver capabilities in itself but hands the PDE over to 
the PDE solver library defined through the \Domain of the PDE.
The general interface is provided through the \LinearPDE class. The
\AdvectivePDE which is derived from the \LinearPDE class
provides an interface to PDE dominated by its advective terms. The \Poisson
class which is also derived form the \LinearPDE class should be used
to define the Poisson equation \index{Poisson}.  

\section{\LinearPDE Class}
\label{SEC LinearPDE}

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:
\begin{equation}\label{LINEARPDE.SINGLE.1}
-(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 \; .
\end{equation}
$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 
\Function on the PDE or objects that can be converted into such \Data objects. 
$A$ is a \RankTwo, $B$, $C$ and $X$ are \RankOne and $D$ and $Y$ are scalar. 
The following natural
boundary conditions are considered \index{boundary condition!natural} on $\Gamma$:
\begin{equation}\label{LINEARPDE.SINGLE.2}
n\hackscore{j}(A\hackscore{jl} u\hackscore{,l}+B\hackscore{j} u)+d u=n\hackscore{j}X\hackscore{j} + y  \;.
\end{equation}
Notice that the coefficients $A$, $B$ and $X$ are defined in the PDE. The coefficients $d$ and $y$ are  
each a \Scalar in the \FunctionOnBoundary.  Constraints \index{constraint} for the solution prescribing the value of the 
solution at certain locations in the domain. They have the form
\begin{equation}\label{LINEARPDE.SINGLE.3}
u=r \mbox{ where } q>0
\end{equation}
$r$ and $q$ are each \Scalar where $q$ is the characteristic function
\index{characteristic function} defining where the constraint is applied.
The constraints defined by \eqn{LINEARPDE.SINGLE.3} override any other condition set by \eqn{LINEARPDE.SINGLE.1}
or \eqn{LINEARPDE.SINGLE.2}. 

For a system of PDEs and a solution with several components the PDE has the form
\begin{equation}\label{LINEARPDE.SYSTEM.1}
-(A\hackscore{ijkl} u\hackscore{k,l}){,j}+(B\hackscore{ijk} u\hackscore{k})\hackscore{,j}+C\hackscore{ikl} u\hackscore{k,l}+D\hackscore{ik} u\hackscore{k} =-X\hackscore{ij,j}+Y\hackscore{i} \; .
\end{equation}
$A$ is a \RankFour, $B$ and $C$ are each a \RankThree, $D$ and $X$ are each a \RankTwo and $Y$ is a \RankOne. 
The natural boundary conditions \index{boundary condition!natural} take the form:
\begin{equation}\label{LINEARPDE.SYSTEM.2}
n\hackscore{j}(A\hackscore{ijkl} u\hackscore{k,l}+B\hackscore{ijk} u\hackscore{k})+d\hackscore{ik} u\hackscore{k}=n\hackscore{j}X\hackscore{ij}+y\hackscore{i}  \;.
\end{equation}
The coefficient $d$ is a \RankTwo and $y$ is a  
\RankOne both in the \FunctionOnBoundary. Constraints \index{constraint} take the form
\begin{equation}\label{LINEARPDE.SYSTEM.3}
u\hackscore{i}=r\hackscore{i} \mbox{ where } q\hackscore{i}>0
\end{equation}
$r$ and $q$ are each \RankOne. Notice that not necessarily all components must 
have a constraint at all locations.

\LinearPDE also supports solution discontinuities \index{discontinuity} 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 
\begin{equation}\label{LINEARPDE.SYSTEM.5}
J\hackscore{ij}=A\hackscore{ijkl}u\hackscore{k,l}+B\hackscore{ijk}u\hackscore{k}-X\hackscore{ij}
\end{equation}
For the case of single solution component and single PDE $J$ is defined
\begin{equation}\label{LINEARPDE.SINGLE.5}
J\hackscore{j}=A\hackscore{jl}u\hackscore{,l}+B\hackscore{j}u\hackscore{k}-X\hackscore{j}
\end{equation}
In the context of discontinuities \index{discontinuity} $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
\begin{equation}\label{LINEARPDE.SYSTEM.6}
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} \; .
\end{equation}
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 \RankTwo and $y^{contact}$ is a  
\RankOne both in the \FunctionOnContactZero or \FunctionOnContactOne.
In case of a single PDE and a single component solution the contact condition takes the form
\begin{equation}\label{LINEARPDE.SINGLE.6}
n\hackscore{j} J^{0}\hackscore{j}=n\hackscore{j} J^{1}\hackscore{j}=y^{contact} - d^{contact}[u]
\end{equation}
In this case the the coefficient $d^{contact}$ and $y^{contact}$ are eaach \Scalar
both in the \FunctionOnContactZero or \FunctionOnContactOne.

The PDE is symmetrical \index{symmetrical} if
\begin{equation}\label{LINEARPDE.SINGLE.4}
A\hackscore{jl}=A\hackscore{lj} \mbox{ and } B\hackscore{j}=C\hackscore{j}
\end{equation}
The system of PDEs is symmetrical \index{symmetrical} if
\begin{eqnarray}
\label{LINEARPDE.SYSTEM.4}
A\hackscore{ijkl}=A\hackscore{klij} \\
B\hackscore{ijk}=C\hackscore{kij} \\
D\hackscore{ik}=D\hackscore{ki} \\
d\hackscore{ik}=d\hackscore{ki} \\
d^{contact}\hackscore{ik}=d^{contact}\hackscore{ki} 
\end{eqnarray}
Note that different from the scalar case~\eqn{LINEARPDE.SINGLE.4} now the coefficients $D$, $d$ abd $d^{contact}$
have to be inspected.

\section{\LinearPDE class}
This is the general class to define a linear PDE in \escript. We list a selction of the most 
important methods of the class only and refer to reference guide \ReferenceGuide for a complete list.

\begin{classdesc}{LinearPDE}{domain,numEquations=0,numSolutions=0}
opens a linear, steady, second order PDE on the \Domain \var{domain}. \var{numEquations}
and \var{numSolutions} gives the number of equations and the number of solutiopn components.
If \var{numEquations} and \var{numSolutions} is non-positive, the number of equations 
and the number solutions, respctively, stay undefined until a coefficient is
defined. 
\end{classdesc}

\begin{methoddesc}[LinearPDE]{setValue}{
\optional{A=Data()}\optional{, B=Data()},
\optional{, C=Data()}\optional{, D=Data()}
\optional{, X=Data()}\optional{, Y=Data()}
\optional{, d=Data()}\optional{, y=Data()}
\optional{, d_contact=Data()}\optional{, y_contact=Data()}
\optional{, q=Data()}\optional{, r=Data()}}
assigns new values to coefficients. 
If the new coefficient value is not a \Data object, it is converted into a \Data object in the
appropriate \FunctionSpace.
\end{methoddesc}

\begin{methoddesc}[LinearPDE]{getCoefficient}{name}
return the value assigned to coefficient \var{name}. If \var{name} is not a valid name 
an exception is raised. 
\end{methoddesc}

\begin{methoddesc}[LinearPDE]{getShapeOfCoefficient}{name}
returns the shape of coefficient \var{name} even if no value has been assigned to it.
\end{methoddesc}

\begin{methoddesc}[LinearPDE]{getFunctionSpaceForCoefficient}{name}
returns the \FunctionSpace of coefficient \var{name} even if no value has been assigned to it.
\end{methoddesc}

\begin{methoddesc}[LinearPDE]{setDebugOn}{}
switches the debug mode to on.
\end{methoddesc}

\begin{methoddesc}[LinearPDE]{setDebugOff}{}
switches the debug mode to on.
\end{methoddesc}

\begin{methoddesc}[LinearPDE]{isUsingLumping}{}
returns \True if \LUMPING is set as the solver for the system of lienar equations.
Otherwise \False is returned.
\end{methoddesc}

\begin{methoddesc}[LinearPDE]{setSolverMethod}{\optional{solver=LinearPDE.DEFAULT}\optional{, preconditioner=LinearPDE.DEFAULT}}
sets the solver method and preconditioner to be used. It is pointed out that a PDE solver library
may not know the specified solver method but may choose a similar method and preconditioner. 
\end{methoddesc}

\begin{methoddesc}[LinearPDE]{getSolverMethodName}{}
returns the name of the solver method and preconditioner which is currently been used.
\end{methoddesc}

\begin{methoddesc}[LinearPDE]{getSolverMethod}{}
returns the solver method and preconditioner which is currently been used.
\end{methoddesc}

\begin{methoddesc}[LinearPDE]{setSolverPackage}{\optional{package=LinearPDE.DEFAULT}}
Set the solver package to be used by PDE library to solve the linear systems of equations. The 
specified package may not be supported by the PDE solver library. In this case, dependng on
the PDE solver, the default solver is used or an exeption is thrown.
If \var{package} is not specified, the default package of the PDE solver library is used. 
\end{methoddesc}

\begin{methoddesc}[LinearPDE]{getSolverPackage}{}
returns the linear solver package currently by the PDE solver library
\end{methoddesc}


\begin{methoddesc}[LinearPDE]{setTolerance}{\optional{tol=1.e-8}}:
resets the tolerance for solution. The actually meaning of tolerance is
depending on the underlying PDE library. In most cases, the tolerance 
will only consider the error from solving the discerete problem but will
not consider any discretization error.
\end{methoddesc}

\begin{methoddesc}[LinearPDE]{getTolerance}{}
returns the current tolerance of the solution
\end{methoddesc}

\begin{methoddesc}[LinearPDE]{getDomain}{}
returns the \Domain of the PDE.
\end{methoddesc}

\begin{methoddesc}[LinearPDE]{getDim}{}
returns the spatial dimension of the PDE.
\end{methoddesc}

\begin{methoddesc}[LinearPDE]{getNumEquations}{}
returns the number of equations.
\end{methoddesc}

\begin{methoddesc}[LinearPDE]{getNumSolutions}{}
returns the number of components of the solution.
\end{methoddesc}

\begin{methoddesc}[LinearPDE]{checkSymmetry}{verbose=\False}
returns \True if the PDE is symmetric and \False otherwise. 
The method is very computational expensive and should only be 
called for testing purposes. The symmetry flag is not altered.
If \var{verbose}=\True information about where symmetry is violated
are printed.
\end{methoddesc}

\begin{methoddesc}[LinearPDE]{getFlux}{u}
returns the flux $J\hackscore{ij}$ \index{flux} for given solution \var{u}
defined by \eqn{LINEARPDE.SYSTEM.5} and \eqn{LINEARPDE.SINGLE.5}, respectively.
\end{methoddesc}


\begin{methoddesc}[LinearPDE]{isSymmetric}{}
returns \True if the PDE has been indicated to be symmetric.
Otherwise \False is returned.
\end{methoddesc}

\begin{methoddesc}[LinearPDE]{setSymmetryOn}{}
indicates that the PDE is symmetric.
\end{methoddesc}

\begin{methoddesc}[LinearPDE]{setSymmetryOff}{}
indicates that the PDE is not symmetric.
\end{methoddesc}

\begin{methoddesc}[LinearPDE]{setReducedOrderOn}{}
switches on the reduction of polynomial order for the solution and equation evaluation even if 
a quadratic or higher interpolation order is defined in the \Domain. This feature may not 
be supported by all PDE libraries.
\end{methoddesc}

\begin{methoddesc}[LinearPDE]{setReducedOrderOff}{}
switches off the reduction of polynomial order for the solution and 
equation evaluation.
\end{methoddesc}

\begin{methoddesc}[LinearPDE]{getOperator}{}
returns the \Operator of the PDE.
\end{methoddesc}

\begin{methoddesc}[LinearPDE]{getRightHandSide}{}
returns the right hand side of the PDE as a \Data object. If
\var{ignoreConstraint}=\True the constraints are not considered
when building up the right hand side.
\end{methoddesc}

\begin{methoddesc}[LinearPDE]{getSystem}{}
returns the \Operator and right hand side of the PDE.
\end{methoddesc}

\begin{methoddesc}[LinearPDE]{getSolution}{
\optional{verbose=False}
\optional{, reordering=LinearPDE.NO_REORDERING}
\optional{, iter_max=1000}
\optional{, drop_tolerance=0.01}
\optional{, drop_storage=1.20}
\optional{, truncation=-1}
\optional{, restart=-1}
}
returns (an approximation of) the solution of the PDE. If \code{verbose=\True} some information during the solution process printed.
\var{reordering} selects a reordering methods that is applied before or during the solution process 
(=\NOREORDERING ,\MINIMUMFILLIN ,\NESTEDDESCTION).
\var{iter_max} specifies the maximum number of iteration steps that are allowed to reach the specified tolerance. 
\var{drop_tolerance} specifies a relative tolerance for small elements to be dropped when building a preconditioner
(eg. in \ILUT). \var{drop_storage} limits the extra storage allowed when building a preconditioner
(eg. in \ILUT). The extra storage is given relative to the size of the stiffness matrix, eg. 
\var{drop_storage=1.2} will allow the preconditioner to use the $1.2$ fold storage space than used 
for the stiffness matrix. \var{truncation} defines the truncation. 
\end{methoddesc}

\begin{memberdesc}[LinearPDE]{DEFAULT}
default method, preconditioner or package to be used to solve the PDE. An appropriate method should be 
chosen by the used PDE solver library.
\end{memberdesc}

\begin{memberdesc}[LinearPDE]{SCSL}
the SCSL library by SGI,~\Ref{SCSL}\footnote{The SCSL library will only be available on SGI systems} 
\end{memberdesc}

\begin{memberdesc}[LinearPDE]{MKL}
the MKL library by Intel,~\Ref{MKL}\footnote{The MKL library will only be available when the intel compilation environment is used.}.
\end{memberdesc}

\begin{memberdesc}[LinearPDE]{UMFPACK}
the UMFPACK,~\Ref{UMFPACK}. Remark: UMFPACK is not parallelized.
\end{memberdesc}

\begin{memberdesc}[LinearPDE]{PASO}
the solver library of \finley, see \Sec{CHAPTER ON FINLEY}.
\end{memberdesc}

\begin{memberdesc}[LinearPDE]{ITERATIVE}
the default iterative method and preconditioner. The actually used method depends on the
PDE solver library and the solver package been choosen. Typically, \PCG is used for symmetric PDEs 
and \BiCGStab otherwise, both with \JACOBI preconditioner. 
\end{memberdesc}

\begin{memberdesc}[LinearPDE]{DIRECT}
the default direct linear solver. 
\end{memberdesc}

\begin{memberdesc}[LinearPDE]{CHOLEVSKY}
direct solver based on Cholevsky factorization (or similar), see~\Ref{Saad}. The solver will require a symmetric PDE. 
\end{memberdesc}

\begin{memberdesc}[LinearPDE]{PCG}
preconditioned conjugate gradient method, see~\Ref{WEISS}\index{linear solver!PCG}\index{PCG}. The solver will require a symmetric PDE.
\end{memberdesc}

\begin{memberdesc}[LinearPDE]{GMRES}
the GMRES method, see~\Ref{WEISS}\index{linear solver!GMRES}\index{GMRES}. Truncation and restart are controlled by the parameters
\var{truncation} and \var{restart} of \method{getSolution}.
\end{memberdesc}

\begin{memberdesc}[LinearPDE]{LUMPING}
uses lumping to solve the system of linear equations~\index{linear solver!lumping}\index{lumping}. This solver technique 
condenses the stiffness matrix to a diagonal matrix so the solution of the linear systems becomes very cheap. It can be used when 
only \var{D} is present but in any case has to applied with care. The difference in the solutions with and without lumping can be significant
but is expect to converge to zero when the mesh gets finer.  
Lumping does not use the linear system solver library.   
\end{memberdesc}

\begin{memberdesc}[LinearPDE]{PRES20}
the GMRES method with truncation after five residuals and
restart after 20 steps, see~\Ref{WEISS}.
\end{memberdesc}[LinearPDE]{CR}

\begin{memberdesc}[LinearPDE]{CGS}
conjugate gradient squared method, see~\Ref{WEISS}.
\end{memberdesc}

\begin{memberdesc}[LinearPDE]{BICGSTAB}
stabilized bi-conjugate gradients methods, see~\Ref{WEISS}. 
\end{memberdesc}

\begin{memberdesc}[LinearPDE]{SSOR}
symmetric successive over-relaxation method, see~\Ref{WEISS}. Typically used as preconditioner but some linear solver libraries support
this as a solver.  
\end{memberdesc}
\begin{memberdesc}[LinearPDE]{ILU0}
the incomplete LU factorization preconditioner with no fill-in, see~\Ref{Saad}. 
\end{memberdesc}

\begin{memberdesc}[LinearPDE]{ILUT}
the incomplete LU factorization preconditioner with fill-in, see~\Ref{Saad}. During the  LU-factorization element with
relative size less then \var{drop_tolerance} are dropped. Moreover, the size of the LU-factorization is restricted to the
\var{drop_storage}-fold of the stiffness matrix. \var{drop_tolerance} and \var{drop_storage} are both set in the 
\method{getSolution} call. 
\end{memberdesc}

\begin{memberdesc}[LinearPDE]{JACOBI}
the Jacobi preconditioner, see~\Ref{Saad}.
\end{memberdesc}

\begin{memberdesc}[LinearPDE]{AMG}
the algebraic--multi grid method, see~\Ref{AMG}. This method can be used as linear solver method but is more robust when used 
in a preconditioner.
\end{memberdesc}

\begin{memberdesc}[LinearPDE]{RILU}
recursive incomplete LU factorization preconditioner, see~\Ref{RILU}. This method is similar to \ILUT but uses smoothing
between levels. During the  LU-factorization element with
relative size less then \var{drop_tolerance} are dropped. Moreover, the size of the LU-factorization is restricted to the
\var{drop_storage}-fold of the stiffness matrix. \var{drop_tolerance} and \var{drop_storage} are both set in the 
\method{getSolution} call. 
\end{memberdesc}

\begin{memberdesc}[LinearPDE]{NO_REORDERING}
no ordering is used during factorization.
\end{memberdesc}

\begin{memberdesc}[LinearPDE]{MINIMUM_FILL_IN}
applies reordering before factorization using a fill-in minimization strategy. You have to check with the particular solver library or
linear solver package if this is supported. In any case, it is advisable to apply reordering on the mesh to minimize fill-in.
\end{memberdesc}

\begin{memberdesc}[LinearPDE]{NESTED_DISSECTION}
applies reordering before factorization using a nested dissection strategy. You have to check with the particular solver library or
linear solver package if this is supported. In any case, it is advisable to apply reordering on the mesh to minimize fill-in.
\end{memberdesc}

\section{The \Poisson Class}
The \Poisson class provides an easy way to define and solve the Poisson
equation
\begin{equation}\label{POISSON.1}
-u\hackscore{,ii}=f\; .
\end{equation}
with homogeneous boundary conditions
\begin{equation}\label{POISSON.2}
n\hackscore{i}u\hackscore{,i}=0
\end{equation}
and homogeneous constraints
\begin{equation}\label{POISSON.3}
u=0 \mbox{ where } q>0
\end{equation}
$f$ has to be a \Scalar in the \Function and $q$ must be
a \Scalar in  the \SolutionFS. 

\begin{classdesc}{Poisson}{domain}
opens a Poisson equation on the \Domain domain. \Poisson is derived from \LinearPDE.
\end{classdesc}
\begin{methoddesc}[Poisson]{setValue}{f=escript.Data(),q=escript.Data()}
assigns new values to \var{f} and \var{q}.
\end{methoddesc}

\section{The \Helmholtz Class}

\section{The \Lame Class}

Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision