1 
jgs 
102 
% $Id$ 
2 



3 


\chapter{The module \linearPDEsPack} 
4 



5 


\declaremodule{extension}{linearPDEs} \modulesynopsis{Linear partial pifferential equation handler} 
6 


The module \linearPDEsPack provides an interface to define and solve linear partial 
7 


differential equations within \escript. \linearPDEsPack does not provide any 
8 


solver capabilities in itself but hand the PDE over to 
9 


the PDE solver library defined through the \Domain of the PDE. 
10 


The general interface is provided through the \LinearPDE class. The 
11 


\AdvectivePDE which is derived from the \LinearPDE class 
12 


provides an interface to PDE dominated by its advective terms. The \Poisson 
13 


class which is also derived form the \LinearPDE class should be used 
14 


to define the Poisson equation \index{Poisson}. 
15 



16 


\section{\LinearPDE Class} 
17 


\label{SEC LinearPDE} 
18 



19 


The \LinearPDE class is used to define a general linear, steady, second order PDE 
20 


for an unknown function $u$ on a given $\Omega$ defined through a \Domain object. 
21 


In the following $\Gamma$ denotes the boundary of the domain $\Omega$. $n$ denotes 
22 


the outer normal field on $\Gamma$. 
23 



24 


A single PDE with a solution with a single component the linear PDE is defined in the 
25 


following form: 
26 


\begin{equation}\label{LINEARPDE.SINGLE.1} 
27 


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


\end{equation} 
29 


$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. 
30 


The coefficients $A$, $B$, $C$, $D$, $X$ and $Y$ have to be specified through \Data objects in the 
31 


\Function on the PDE or objects that can be converted into such \Data objects. 
32 


$A$ is a \RankTwo, $B$, $C$ and $X$ are \RankOne and $D$ and $Y$ are scalar. 
33 


The following natural 
34 


boundary conditions are considered \index{boundary condition!natural} on $\Gamma$: 
35 


\begin{equation}\label{LINEARPDE.SINGLE.2} 
36 


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


\end{equation} 
38 


Notice that the coefficients $A$, $B$ and $X$ are defined in the PDE. The coefficients $d$ and $y$ are 
39 


each a \Scalar in the \FunctionOnBoundary. Constraints \index{constraint} for the solution prescribing the value of the 
40 


solution at certain locations in the domain. They have the form 
41 


\begin{equation}\label{LINEARPDE.SINGLE.3} 
42 


u=r \mbox{ where } q>0 
43 


\end{equation} 
44 


$r$ and $q$ are each \Scalar where $q$ is the characteristic function 
45 


\index{characteristic function} defining where the constraint is applied. 
46 


The constraints defined by \eqn{LINEARPDE.SINGLE.3} override any other condition set by \eqn{LINEARPDE.SINGLE.1} 
47 


or \eqn{LINEARPDE.SINGLE.2}. The PDE is symmetrical \index{symmetrical} if 
48 


\begin{equation}\label{LINEARPDE.SINGLE.4} 
49 


A\hackscore{jl}=A\hackscore{lj} \mbox{ and } B\hackscore{j}=C\hackscore{j} 
50 


\end{equation} 
51 


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


\begin{equation}\label{LINEARPDE.SYSTEM.1} 
53 


(A\hackscore{ijkl} u\hackscore{k,l}){,j}+(B\hackscore{ijk} u_k)\hackscore{,j}+C\hackscore{ikl} u\hackscore{k,l}+D\hackscore{ik} u_k =X\hackscore{ij,j}+Y\hackscore{i} \; . 
54 


\end{equation} 
55 


$A$ is a \RankFour, $B$ and $C$ are each a \RankThree, $D$ and $X$ are each a \RankTwo and $Y$ is a \RankOne. 
56 


The natural boundary conditions \index{boundary condition!natural} take the form: 
57 


\begin{equation}\label{LINEARPDE.SYSTEM.2} 
58 


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


\end{equation} 
60 


The coefficient $d$ is a \RankTwo and $y$ is a 
61 


\RankOne both in the \FunctionOnBoundary. Constraints \index{constraint} take the form 
62 


\begin{equation}\label{LINEARPDE.SYSTEM.3} 
63 


u\hackscore{i}=r\hackscore{i} \mbox{ where } q\hackscore{i}>0 
64 


\end{equation} 
65 


$r$ and $q$ are each \RankOne. Notice that at some locations not necessarily all components must 
66 


have a constraint. The system of PDEs is symmetrical \index{symmetrical} if 
67 


\begin{eqnarray}\label{LINEARPDE.SYSTEM.4} 
68 


A\hackscore{ijkl}=A\hackscore{klij} \\ 
69 


B\hackscore{ijk}=C\hackscore{kij} \\ 
70 


D\hackscore{ik}=D\hackscore{ki} \\ 
71 


d\hackscore{ik}=d\hackscore{ki} \ 
72 


\end{eqnarray} 
73 


\LinearPDE also supports solution discontinuities \index{discontinuity} over contact region $\Gamma^{contact}$ 
74 


in the domain $\Omega$. To specify the conditions across the discontinuity we are using the 
75 


generalised flux $J$ which is in the case of a systems of PDEs and several components of the solution 
76 


defined as 
77 


\begin{equation}\label{LINEARPDE.SYSTEM.5} 
78 


J\hackscore{ij}=A\hackscore{ijkl}u\hackscore{k,l}+B\hackscore{ijk}u\hackscore{k}X\hackscore{ij} 
79 


\end{equation} 
80 


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


\begin{equation}\label{LINEARPDE.SINGLE.5} 
82 


J\hackscore{j}=A\hackscore{jl}u\hackscore{,l}+B\hackscore{j}u\hackscore{k}X\hackscore{j} 
83 


\end{equation} 
84 


In the context of discontinuities \index{discontinuity} $n$ denotes the normal on the 
85 


discontinuity pointing from side 0 towards side 1. For a system of PDEs 
86 


the contact condition takes the form 
87 


\begin{equation}\label{LINEARPDE.SYSTEM.6} 
88 


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} \; . 
89 


\end{equation} 
90 


where $J^{0}$ and $J^{1}$ are the fluxes on side $0$ and side $1$ of the 
91 


discontinuity $\Gamma^{contact}$, respectively. $[u]$, which is the difference 
92 


of the solution at side 1 and at side 0, denotes the jump of $u$ across $\Gamma^{contact}$. 
93 


The coefficient $d^{contact}$ is a \RankTwo and $y^{contact}$ is a 
94 


\RankOne both in the \FunctionOnContactZero or \FunctionOnContactOne. 
95 


In case of a single PDE and a single component solution the contact condition takes the form 
96 


\begin{equation}\label{LINEARPDE.SINGLE.6} 
97 


n\hackscore{j} J^{0}\hackscore{j}=n\hackscore{j} J^{1}\hackscore{j}=y^{contact}  d^{contact}[u] 
98 


\end{equation} 
99 


In this case the the coefficient $d^{contact}$ and $y^{contact}$ are eaach \Scalar 
100 


both in the \FunctionOnContactZero or \FunctionOnContactOne. 
101 



102 


\begin{classdesc}{LinearPDE}{domain,numEquations=0,numSolutions=0} 
103 


opens a linear, steady, second order PDE on the \Domain \var{domain}. \var{numEquations} 
104 


and \var{numSolutions} gives the number of equations and the number of solutiopn components. 
105 


If \var{numEquations} and \var{numSolutions} is nonpositive, the number of equations 
106 


and the number solutions, respctively, stay undefined until a coefficient is 
107 


defined. 
108 


\end{classdesc} 
109 



110 


\begin{methoddesc}[LinearPDE]{setValues}{arg1,arg2,...,argN} 
111 


assigns new values to coefficients \var{arg1}, \var{arg2}, $\ldots$ \var{argN}. 
112 


The coefficients must be one of the valid coefficient names 
113 


\var{A}, 
114 


\var{B}, 
115 


\var{C}, 
116 


\var{D}, 
117 


\var{X}, 
118 


\var{Y}, 
119 


\var{d}, 
120 


\var{y}, 
121 


\var{dcontact}, 
122 


\var{ycontact}, 
123 


\var{r}, 
124 


or \var{q}. 
125 


If the coefficient is not a \Data object, it is converted into a \Data object in the 
126 


appropriate \FunctionSpace. 
127 


\end{methoddesc} 
128 



129 


\begin{methoddesc}[LinearPDE]{getCoefficient}{name} 
130 


return the value assigned to coefficient \var{name}. If \var{name} is not a valid name 
131 


an exception is raised. 
132 


\end{methoddesc} 
133 



134 


\begin{methoddesc}[LinearPDE]{cleanCoefficients}{} 
135 


resets all coefficients to their initialization values. This method is useful to call when trying to save memory 
136 


as is releaces eqnerences to coefficients but keeping the differential operator. 
137 


\end{methoddesc} 
138 



139 


\begin{methoddesc}[LinearPDE]{getShapeOfCoefficient}{name} 
140 


returns the shape of coefficient \var{name} even if no value has been assigned to it. 
141 


\end{methoddesc} 
142 



143 



144 


\begin{methoddesc}[LinearPDE]{getFunctionSpaceOfCoefficient}{name} 
145 


returns the \FunctionSpace of coefficient \var{name} even if no value has been assigned to it. 
146 


\end{methoddesc} 
147 



148 


\begin{methoddesc}[LinearPDE]{hasCoefficient}{name} 
149 


returns \True if \var{name} is valid name of a coefficient 
150 


\end{methoddesc} 
151 



152 



153 


\begin{methoddesc}[LinearPDE]{getFunctionSpaceForEquation}{} 
154 


returns \FunctionSpace of the right hand side which is identical to \FunctionSpace of 
155 


the result of the PDE operator. 
156 


\end{methoddesc} 
157 



158 


\begin{methoddesc}[LinearPDE]{getFunctionSpaceForSolution}{} 
159 


returns \FunctionSpace of the solution of the PDE hich is identical to \FunctionSpace of 
160 


the result of argument of the PDE operator. 
161 


\end{methoddesc} 
162 



163 


\begin{methoddesc}[LinearPDE]{setDebugOn}{} 
164 


switches the debug mode to on. 
165 


\end{methoddesc} 
166 



167 


\begin{methoddesc}[LinearPDE]{setDebugOff}{} 
168 


switches the debug mode to on. 
169 


\end{methoddesc} 
170 



171 


\begin{methoddesc}[LinearPDE]{debug}{} 
172 


returns \True if the debug mode is switched on. Otherwise it return \False. 
173 


\end{methoddesc} 
174 



175 


\begin{methoddesc}[LinearPDE]{setLumpingOn}{} 
176 


switches on lumping \index{lumping}. If lumping is switched on 
177 


the operator is 
178 


\end{methoddesc} 
179 



180 



181 


\begin{methoddesc}[LinearPDE]{setLumpingOff} 
182 


switches lumping off \index{lumping}. 
183 


\end{methoddesc} 
184 



185 


\begin{methoddesc}[LinearPDE]{setLumping}{flag=\False} 
186 


switches on lumping if \var{flag} is \True. Otherwise lumping is swiched off. 
187 


\end{methoddesc} 
188 



189 


\begin{methoddesc}[LinearPDE]{isUsingLumping}{} 
190 


returns \True if lumping is switched on. Otherwise \False is returned. 
191 


\end{methoddesc} 
192 



193 


\begin{memberdesc}[LinearPDE]{DEFAULT_METHOD} 
194 


default method to be used to solve the PDE. An appropriate method should be 
195 


chosen by the used PDE solver library. 
196 


\end{memberdesc} 
197 



198 


\begin{memberdesc}[LinearPDE]{DIRECT} 
199 


direct linear solver~\Ref{SAAD} 
200 


\end{memberdesc} 
201 



202 


\begin{memberdesc}[LinearPDE]{CHOLEVSKY} 
203 


direct solver based on Cholevsky factorization (or similar), see~\Ref{SAAD}. The solver will require a symmetric PDE. 
204 


\end{memberdesc} 
205 



206 


\begin{memberdesc}[LinearPDE]{PCG} 
207 


preconditioned conjugate gradient method, see~\Ref{WEISS}. The solver will require a symmetric PDE. 
208 


\end{memberdesc} 
209 



210 


\begin{memberdesc}[LinearPDE]{GMRES} 
211 


the GMRES method, see~\Ref{WEISS}. Truncation and restart ar econtrolled by the parameters 
212 


\var{truncation} and \var{restart} of \method{getSolution}. 
213 


\end{memberdesc} 
214 



215 


\begin{memberdesc}[LinearPDE]{PRES20} 
216 


the GMRES method with trunction after five residuals and 
217 


restart after 20 steps, see~\Ref{WEISS}. 
218 


\end{memberdesc} 
219 



220 


\begin{memberdesc}[LinearPDE]{CR} 
221 


conjugate residual method, see~\Ref{WEISS}. 
222 


\end{memberdesc} 
223 



224 


\begin{memberdesc}[LinearPDE]{CGS} 
225 


conjugate gradient squared method, see~\Ref{WEISS}. 
226 


\end{memberdesc} 
227 



228 


\begin{memberdesc}[LinearPDE]{BICGSTAB} 
229 


stabilzed biconjugate gradients methods, see~\Ref{WEISS}. 
230 


\end{memberdesc} 
231 



232 


\begin{memberdesc}[LinearPDE]{SSOR} 
233 


symmetric successive overrelaxtion method, see~\Ref{WEISS}. 
234 


\end{memberdesc} 
235 



236 


\begin{methoddesc}[LinearPDE]{setSolverMethod}{solver=linearPDE.DEFAULT_METHOD} 
237 


sets the solver method to be used. It is pointed out that the PDE solver library 
238 


does not know the specified solver method but may choose a similar method. 
239 


\end{methoddesc} 
240 



241 


\begin{methoddesc}[LinearPDE]{setTolerance}{tol=1.e8} 
242 


resets the tolerance for solution. The actually meaning of tolerance is 
243 


depending on the underlying PDE library. In most cases, the tolerance 
244 


will only consider the error from solving the discerete problem but will 
245 


not consider any discretization error. 
246 


\end{methoddesc} 
247 



248 


\begin{methoddesc}[LinearPDE]{getTolerance}{} 
249 


returns the current tolerance of the solution 
250 


\end{methoddesc} 
251 



252 


\begin{methoddesc}[LinearPDE]{isSymmetric}{} 
253 


returns \True if the PDE has been indicated to be symmetric. 
254 


Otherwise \False is returned. 
255 


\end{methoddesc} 
256 



257 


\begin{methoddesc}[LinearPDE]{setSymmetryOn}{} 
258 


indicates that the PDE is symmetric. 
259 


\end{methoddesc} 
260 



261 


\begin{methoddesc}[LinearPDE]{setSymmetryOff}{} 
262 


indicates that the PDE is not symmetric. 
263 


\end{methoddesc} 
264 



265 


\begin{methoddesc}[LinearPDE]{setSymmetryTo}{flag=\False} 
266 


indicates that the PDE is symmetric if \var{flag}=\True 
267 


and indicates a nonsymmetric PDE is \var{flag}=\False. 
268 


\end{methoddesc} 
269 



270 


\begin{methoddesc}[LinearPDE]{setReducedOrderOn}{} 
271 


switches on the reduction of polynomial order for the solution and 
272 


equation evaluation. 
273 


\end{methoddesc} 
274 



275 


\begin{methoddesc}[LinearPDE]{setReducedOrderOff}{} 
276 


switches off the reduction of polynomial order for the solution and 
277 


equation evaluation. 
278 


\end{methoddesc} 
279 



280 


\begin{methoddesc}[LinearPDE]{setReducedOrderTo}{flag=\False} 
281 


switches on the reduction of polynomial order for the solution and 
282 


equation evaluation if \var{flag}=\True. Otherwise 
283 


the order reduction is switched off. 
284 


\end{methoddesc} 
285 



286 


\begin{methoddesc}[LinearPDE]{setReducedOrderForSolutionOn}{} 
287 


switches on reduction of polynomial order for the solution. 
288 


\end{methoddesc} 
289 



290 


\begin{methoddesc}[LinearPDE]{setReducedOrderForSolutionOff}{} 
291 


switches off reduction of polynomial order for the solution. 
292 


\end{methoddesc} 
293 



294 


\begin{methoddesc}[LinearPDE]{setReducedOrderForSolutionTo}{flag=\False} 
295 


switches on the reduction of polynomial order for the solution 
296 


if \var{flag}=\True. Otherwise 
297 


the order reduction is switched off. 
298 


\end{methoddesc} 
299 



300 


\begin{methoddesc}[LinearPDE]{setReducedOrderForEquationOn}{} 
301 


switches on reduction of polynomial order for the equation evaluation. 
302 


\end{methoddesc} 
303 



304 


\begin{methoddesc}[LinearPDE]{setReducedOrderForEquationOff}{} 
305 


switches off reduction of polynomial order for the equation evaluation. 
306 


\end{methoddesc} 
307 



308 


\begin{methoddesc}[LinearPDE]{setReducedOrderForEquationTo}{flag=\False} 
309 


switches on the reduction of polynomial order for the equation 
310 


evaluation if \var{flag}=\True. Otherwise 
311 


the order reduction is switched off. 
312 


\end{methoddesc} 
313 



314 


\begin{methoddesc}[LinearPDE]{getOperator}{} 
315 


returns the \Operator of the PDE. 
316 


\end{methoddesc} 
317 



318 


\begin{methoddesc}[LinearPDE]{getRightHandSide}{ignoreConstraint=\False} 
319 


returns the right hand side of the PDE as a \Data object. If 
320 


\var{ignoreConstraint}=\True the constraints are not considered 
321 


when building up the right hand side. 
322 


\end{methoddesc} 
323 



324 


\begin{methoddesc}[LinearPDE]{getSystem}{} 
325 


returns the \Operator and right hand side of the PDE. 
326 


\end{methoddesc} 
327 



328 


\begin{methoddesc}[LinearPDE]{getSolution}{option1,option2,...,optionN} 
329 


returns (an approximation of) the solution of the PDE. \var{options1}, \var{options2} 
330 


$\ldots$ \var{optionsN} are options handed over to the underlying PDE solver library. 
331 


\end{methoddesc} 
332 



333 


\begin{methoddesc}[LinearPDE]{getDomain}{} 
334 


returns the \Domain of the PDE. 
335 


\end{methoddesc} 
336 



337 


\begin{methoddesc}[LinearPDE]{getDim}{} 
338 


returns the spatial dimension of the PDE. 
339 


\end{methoddesc} 
340 



341 


\begin{methoddesc}[LinearPDE]{getNumEquations}{} 
342 


returns the number of equations. 
343 


\end{methoddesc} 
344 



345 


\begin{methoddesc}[LinearPDE]{getNumSolutions}{} 
346 


returns the number of components of the solution. 
347 


\end{methoddesc} 
348 



349 


\begin{methoddesc}[LinearPDE]{checkSymmetry}{verbose=\False} 
350 


returns \True if the PDE is symmetric and \False otherwise. 
351 


The method is very computational expensive and should only be 
352 


called for testing purposes. The symmetry flag is not altered. 
353 


If \var{verbose}=\True information about where symmetry is violated 
354 


are printed. 
355 


\end{methoddesc} 
356 



357 


\begin{methoddesc}[LinearPDE]{getFlux}{u} 
358 


returns the flux $J\hackscore{ij}$ \index{flux} for given solution \var{u} 
359 


defined by \eqn{LINEARPDE.SYSTEM.5} and \eqn{LINEARPDE.SINGLE.5}, respectively. 
360 


\end{methoddesc} 
361 



362 


\begin{methoddesc}[LinearPDE]{applyOperator}{u} 
363 


applies the PDE operator to \var{u} 
364 


\end{methoddesc} 
365 



366 


\begin{methoddesc}[LinearPDE]{getResidual}{u} 
367 


returns the residual when insering \var{u} into the PDE 
368 


\end{methoddesc} 
369 



370 


\section{\AdvectivePDE Class} 
371 


under construction 
372 



373 


\section{The \Poisson Class} 
374 



375 


The \Poisson class provides an easy way to define and solve the Poisson 
376 


equation 
377 


\begin{equation}\label{POISSON.1} 
378 


u\hackscore{,ii}=f\; . 
379 


\end{equation} 
380 


with homogeneous boundary conditions 
381 


\begin{equation}\label{POISSON.2} 
382 


n\hackscore{i}u\hackscore{,i}=0 
383 


\end{equation} 
384 


and homogeneous constraints 
385 


\begin{equation}\label{POISSON.3} 
386 


u=0 \mbox{ where } q>0 
387 


\end{equation} 
388 


$f$ has to be a \Scalar in the \Function and $q$ must be 
389 


a \Scalar in the \SolutionFS. 
390 



391 


\begin{classdesc}{Poisson}{domain} 
392 


opens a Poisson equation on the \Domain domain. \Poisson is derived from \LinearPDE. 
393 


\end{classdesc} 
394 


\begin{methoddesc}[Poisson]{setValue}{f=escript.Data(),q=escript.Data()} 
395 


assigns new values to \var{f} and \var{q}. 
396 


\end{methoddesc} 