 # Contents of /trunk/doc/user/linearPDE.tex

Revision 653 - (show annotations)
Fri Mar 24 07:34:34 2006 UTC (14 years, 8 months ago) by gross
File MIME type: application/x-tex
File size: 19678 byte(s)
linearPDE docu fixed

 1 % $Id$ 2 % 3 % Copyright © 2006 by ACcESS MNRF 4 % \url{http://www.access.edu.au 5 % Primary Business: Queensland, Australia. 6 % Licensed under the Open Software License version 3.0 7 8 % 9 10 11 \chapter{The module \linearPDEs} 12 13 \declaremodule{extension}{linearPDEs} \modulesynopsis{Linear partial pifferential equation handler} 14 The module \linearPDEs provides an interface to define and solve linear partial 15 differential equations within \escript. \linearPDEs does not provide any 16 solver capabilities in itself but hands the PDE over to 17 the PDE solver library defined through the \Domain of the PDE. 18 The general interface is provided through the \LinearPDE class. The 19 \AdvectivePDE which is derived from the \LinearPDE class 20 provides an interface to PDE dominated by its advective terms. The \Poisson 21 class which is also derived form the \LinearPDE class should be used 22 to define the Poisson equation \index{Poisson}. 23 24 \section{\LinearPDE Class} 25 \label{SEC LinearPDE} 26 27 The \LinearPDE class is used to define a general linear, steady, second order PDE 28 for an unknown function $u$ on a given $\Omega$ defined through a \Domain object. 29 In the following $\Gamma$ denotes the boundary of the domain $\Omega$. $n$ denotes 30 the outer normal field on $\Gamma$. 31 32 For a single PDE with a solution with a single component the linear PDE is defined in the 33 following form: 34 \begin{equation}\label{LINEARPDE.SINGLE.1} 35 -(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 \; . 36 \end{equation} 37 $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. 38 The coefficients $A$, $B$, $C$, $D$, $X$ and $Y$ have to be specified through \Data objects in the 39 \Function on the PDE or objects that can be converted into such \Data objects. 40 $A$ is a \RankTwo, $B$, $C$ and $X$ are \RankOne and $D$ and $Y$ are scalar. 41 The following natural 42 boundary conditions are considered \index{boundary condition!natural} on $\Gamma$: 43 \begin{equation}\label{LINEARPDE.SINGLE.2} 44 n\hackscore{j}(A\hackscore{jl} u\hackscore{,l}+B\hackscore{j} u)+d u=n\hackscore{j}X\hackscore{j} + y \;. 45 \end{equation} 46 Notice that the coefficients $A$, $B$ and $X$ are defined in the PDE. The coefficients $d$ and $y$ are 47 each a \Scalar in the \FunctionOnBoundary. Constraints \index{constraint} for the solution prescribing the value of the 48 solution at certain locations in the domain. They have the form 49 \begin{equation}\label{LINEARPDE.SINGLE.3} 50 u=r \mbox{ where } q>0 51 \end{equation} 52 $r$ and $q$ are each \Scalar where $q$ is the characteristic function 53 \index{characteristic function} defining where the constraint is applied. 54 The constraints defined by \eqn{LINEARPDE.SINGLE.3} override any other condition set by \eqn{LINEARPDE.SINGLE.1} 55 or \eqn{LINEARPDE.SINGLE.2}. 56 57 For a system of PDEs and a solution with several components the PDE has the form 58 \begin{equation}\label{LINEARPDE.SYSTEM.1} 59 -(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} \; . 60 \end{equation} 61 $A$ is a \RankFour, $B$ and $C$ are each a \RankThree, $D$ and $X$ are each a \RankTwo and $Y$ is a \RankOne. 62 The natural boundary conditions \index{boundary condition!natural} take the form: 63 \begin{equation}\label{LINEARPDE.SYSTEM.2} 64 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} \;. 65 \end{equation} 66 The coefficient $d$ is a \RankTwo and $y$ is a 67 \RankOne both in the \FunctionOnBoundary. Constraints \index{constraint} take the form 68 \begin{equation}\label{LINEARPDE.SYSTEM.3} 69 u\hackscore{i}=r\hackscore{i} \mbox{ where } q\hackscore{i}>0 70 \end{equation} 71 $r$ and $q$ are each \RankOne. Notice that not necessarily all components must 72 have a constraint at all locations. 73 74 \LinearPDE also supports solution discontinuities \index{discontinuity} over contact region $\Gamma^{contact}$ 75 in the domain $\Omega$. To specify the conditions across the discontinuity we are using the 76 generalised flux $J$ which is in the case of a systems of PDEs and several components of the solution 77 defined as 78 \begin{equation}\label{LINEARPDE.SYSTEM.5} 79 J\hackscore{ij}=A\hackscore{ijkl}u\hackscore{k,l}+B\hackscore{ijk}u\hackscore{k}-X\hackscore{ij} 80 \end{equation} 81 For the case of single solution component and single PDE $J$ is defined 82 \begin{equation}\label{LINEARPDE.SINGLE.5} 83 J\hackscore{j}=A\hackscore{jl}u\hackscore{,l}+B\hackscore{j}u\hackscore{k}-X\hackscore{j} 84 \end{equation} 85 In the context of discontinuities \index{discontinuity} $n$ denotes the normal on the 86 discontinuity pointing from side 0 towards side 1. For a system of PDEs 87 the contact condition takes the form 88 \begin{equation}\label{LINEARPDE.SYSTEM.6} 89 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} \; . 90 \end{equation} 91 where $J^{0}$ and $J^{1}$ are the fluxes on side $0$ and side $1$ of the 92 discontinuity $\Gamma^{contact}$, respectively. $[u]$, which is the difference 93 of the solution at side 1 and at side 0, denotes the jump of $u$ across $\Gamma^{contact}$. 94 The coefficient $d^{contact}$ is a \RankTwo and $y^{contact}$ is a 95 \RankOne both in the \FunctionOnContactZero or \FunctionOnContactOne. 96 In case of a single PDE and a single component solution the contact condition takes the form 97 \begin{equation}\label{LINEARPDE.SINGLE.6} 98 n\hackscore{j} J^{0}\hackscore{j}=n\hackscore{j} J^{1}\hackscore{j}=y^{contact} - d^{contact}[u] 99 \end{equation} 100 In this case the the coefficient $d^{contact}$ and $y^{contact}$ are eaach \Scalar 101 both in the \FunctionOnContactZero or \FunctionOnContactOne. 102 103 The PDE is symmetrical \index{symmetrical} if 104 \begin{equation}\label{LINEARPDE.SINGLE.4} 105 A\hackscore{jl}=A\hackscore{lj} \mbox{ and } B\hackscore{j}=C\hackscore{j} 106 \end{equation} 107 The system of PDEs is symmetrical \index{symmetrical} if 108 \begin{eqnarray} 109 \label{LINEARPDE.SYSTEM.4} 110 A\hackscore{ijkl}=A\hackscore{klij} \\ 111 B\hackscore{ijk}=C\hackscore{kij} \\ 112 D\hackscore{ik}=D\hackscore{ki} \\ 113 d\hackscore{ik}=d\hackscore{ki} \\ 114 d^{contact}\hackscore{ik}=d^{contact}\hackscore{ki} 115 \end{eqnarray} 116 Note that different from the scalar case~\eqn{LINEARPDE.SINGLE.4} now the coefficients $D$, $d$ abd $d^{contact}$ 117 have to be inspected. 118 119 \section{\LinearPDE class} 120 This is the general class to define a linear PDE in \escript. We list a selction of the most 121 important methods of the class only and refer to reference guide \ReferenceGuide for a complete list. 122 123 \begin{classdesc}{LinearPDE}{domain,numEquations=0,numSolutions=0} 124 opens a linear, steady, second order PDE on the \Domain \var{domain}. \var{numEquations} 125 and \var{numSolutions} gives the number of equations and the number of solutiopn components. 126 If \var{numEquations} and \var{numSolutions} is non-positive, the number of equations 127 and the number solutions, respctively, stay undefined until a coefficient is 128 defined. 129 \end{classdesc} 130 131 \begin{methoddesc}[LinearPDE]{setValue}{ 132 \optional{A=Data()}\optional{, B=Data()}, 133 \optional{, C=Data()}\optional{, D=Data()} 134 \optional{, X=Data()}\optional{, Y=Data()} 135 \optional{, d=Data()}\optional{, y=Data()} 136 \optional{, d_contact=Data()}\optional{, y_contact=Data()} 137 \optional{, q=Data()}\optional{, r=Data()}} 138 assigns new values to coefficients. 139 If the new coefficient value is not a \Data object, it is converted into a \Data object in the 140 appropriate \FunctionSpace. 141 \end{methoddesc} 142 143 \begin{methoddesc}[LinearPDE]{getCoefficient}{name} 144 return the value assigned to coefficient \var{name}. If \var{name} is not a valid name 145 an exception is raised. 146 \end{methoddesc} 147 148 \begin{methoddesc}[LinearPDE]{getShapeOfCoefficient}{name} 149 returns the shape of coefficient \var{name} even if no value has been assigned to it. 150 \end{methoddesc} 151 152 \begin{methoddesc}[LinearPDE]{getFunctionSpaceForCoefficient}{name} 153 returns the \FunctionSpace of coefficient \var{name} even if no value has been assigned to it. 154 \end{methoddesc} 155 156 \begin{methoddesc}[LinearPDE]{setDebugOn}{} 157 switches the debug mode to on. 158 \end{methoddesc} 159 160 \begin{methoddesc}[LinearPDE]{setDebugOff}{} 161 switches the debug mode to on. 162 \end{methoddesc} 163 164 \begin{methoddesc}[LinearPDE]{isUsingLumping}{} 165 returns \True if \LUMPING is set as the solver for the system of lienar equations. 166 Otherwise \False is returned. 167 \end{methoddesc} 168 169 \begin{methoddesc}[LinearPDE]{setSolverMethod}{\optional{solver=LinearPDE.DEFAULT}\optional{, preconditioner=LinearPDE.DEFAULT}} 170 sets the solver method and preconditioner to be used. It is pointed out that a PDE solver library 171 may not know the specified solver method but may choose a similar method and preconditioner. 172 \end{methoddesc} 173 174 \begin{methoddesc}[LinearPDE]{getSolverMethodName}{} 175 returns the name of the solver method and preconditioner which is currently been used. 176 \end{methoddesc} 177 178 \begin{methoddesc}[LinearPDE]{getSolverMethod}{} 179 returns the solver method and preconditioner which is currently been used. 180 \end{methoddesc} 181 182 \begin{methoddesc}[LinearPDE]{setSolverPackage}{\optional{package=LinearPDE.DEFAULT}} 183 Set the solver package to be used by PDE library to solve the linear systems of equations. The 184 specified package may not be supported by the PDE solver library. In this case, dependng on 185 the PDE solver, the default solver is used or an exeption is thrown. 186 If \var{package} is not specified, the default package of the PDE solver library is used. 187 \end{methoddesc} 188 189 \begin{methoddesc}[LinearPDE]{getSolverPackage}{} 190 returns the linear solver package currently by the PDE solver library 191 \end{methoddesc} 192 193 194 \begin{methoddesc}[LinearPDE]{setTolerance}{\optional{tol=1.e-8}}: 195 resets the tolerance for solution. The actually meaning of tolerance is 196 depending on the underlying PDE library. In most cases, the tolerance 197 will only consider the error from solving the discerete problem but will 198 not consider any discretization error. 199 \end{methoddesc} 200 201 \begin{methoddesc}[LinearPDE]{getTolerance}{} 202 returns the current tolerance of the solution 203 \end{methoddesc} 204 205 \begin{methoddesc}[LinearPDE]{getDomain}{} 206 returns the \Domain of the PDE. 207 \end{methoddesc} 208 209 \begin{methoddesc}[LinearPDE]{getDim}{} 210 returns the spatial dimension of the PDE. 211 \end{methoddesc} 212 213 \begin{methoddesc}[LinearPDE]{getNumEquations}{} 214 returns the number of equations. 215 \end{methoddesc} 216 217 \begin{methoddesc}[LinearPDE]{getNumSolutions}{} 218 returns the number of components of the solution. 219 \end{methoddesc} 220 221 \begin{methoddesc}[LinearPDE]{checkSymmetry}{verbose=\False} 222 returns \True if the PDE is symmetric and \False otherwise. 223 The method is very computational expensive and should only be 224 called for testing purposes. The symmetry flag is not altered. 225 If \var{verbose}=\True information about where symmetry is violated 226 are printed. 227 \end{methoddesc} 228 229 \begin{methoddesc}[LinearPDE]{getFlux}{u} 230 returns the flux $J\hackscore{ij}$ \index{flux} for given solution \var{u} 231 defined by \eqn{LINEARPDE.SYSTEM.5} and \eqn{LINEARPDE.SINGLE.5}, respectively. 232 \end{methoddesc} 233 234 235 \begin{methoddesc}[LinearPDE]{isSymmetric}{} 236 returns \True if the PDE has been indicated to be symmetric. 237 Otherwise \False is returned. 238 \end{methoddesc} 239 240 \begin{methoddesc}[LinearPDE]{setSymmetryOn}{} 241 indicates that the PDE is symmetric. 242 \end{methoddesc} 243 244 \begin{methoddesc}[LinearPDE]{setSymmetryOff}{} 245 indicates that the PDE is not symmetric. 246 \end{methoddesc} 247 248 \begin{methoddesc}[LinearPDE]{setReducedOrderOn}{} 249 switches on the reduction of polynomial order for the solution and equation evaluation even if 250 a quadratic or higher interpolation order is defined in the \Domain. This feature may not 251 be supported by all PDE libraries. 252 \end{methoddesc} 253 254 \begin{methoddesc}[LinearPDE]{setReducedOrderOff}{} 255 switches off the reduction of polynomial order for the solution and 256 equation evaluation. 257 \end{methoddesc} 258 259 \begin{methoddesc}[LinearPDE]{getOperator}{} 260 returns the \Operator of the PDE. 261 \end{methoddesc} 262 263 \begin{methoddesc}[LinearPDE]{getRightHandSide}{} 264 returns the right hand side of the PDE as a \Data object. If 265 \var{ignoreConstraint}=\True the constraints are not considered 266 when building up the right hand side. 267 \end{methoddesc} 268 269 \begin{methoddesc}[LinearPDE]{getSystem}{} 270 returns the \Operator and right hand side of the PDE. 271 \end{methoddesc} 272 273 \begin{methoddesc}[LinearPDE]{getSolution}{ 274 \optional{verbose=False} 275 \optional{, reordering=LinearPDE.NO_REORDERING} 276 \optional{, iter_max=1000} 277 \optional{, drop_tolerance=0.01} 278 \optional{, drop_storage=1.20} 279 \optional{, truncation=-1} 280 \optional{, restart=-1} 281 } 282 returns (an approximation of) the solution of the PDE. If \code{verbose=\True} some information during the solution process printed. 283 \var{reordering} selects a reordering methods that is applied before or during the solution process 284 (=\NOREORDERING ,\MINIMUMFILLIN ,\NESTEDDESCTION). 285 \var{iter_max} specifies the maximum number of iteration steps that are allowed to reach the specified tolerance. 286 \var{drop_tolerance} specifies a relative tolerance for small elements to be dropped when building a preconditioner 287 (eg. in \ILUT). \var{drop_storage} limits the extra storage allowed when building a preconditioner 288 (eg. in \ILUT). The extra storage is given relative to the size of the stiffness matrix, eg. 289 \var{drop_storage=1.2} will allow the preconditioner to use the $1.2$ fold storage space than used 290 for the stiffness matrix. \var{truncation} defines the truncation. 291 \end{methoddesc} 292 293 \begin{memberdesc}[LinearPDE]{DEFAULT} 294 default method, preconditioner or package to be used to solve the PDE. An appropriate method should be 295 chosen by the used PDE solver library. 296 \end{memberdesc} 297 298 \begin{memberdesc}[LinearPDE]{SCSL} 299 the SCSL library by SGI,~\Ref{SCSL}\footnote{The SCSL library will only be available on SGI systems} 300 \end{memberdesc} 301 302 \begin{memberdesc}[LinearPDE]{MKL} 303 the MKL library by Intel,~\Ref{MKL}\footnote{The MKL library will only be available when the intel compilation environment is used.}. 304 \end{memberdesc} 305 306 \begin{memberdesc}[LinearPDE]{UMFPACK} 307 the UMFPACK,~\Ref{UMFPACK}. Remark: UMFPACK is not parallelized. 308 \end{memberdesc} 309 310 \begin{memberdesc}[LinearPDE]{PASO} 311 the solver library of \finley, see \Sec{CHAPTER ON FINLEY}. 312 \end{memberdesc} 313 314 \begin{memberdesc}[LinearPDE]{ITERATIVE} 315 the default iterative method and preconditioner. The actually used method depends on the 316 PDE solver library and the solver package been choosen. Typically, \PCG is used for symmetric PDEs 317 and \BiCGStab otherwise, both with \JACOBI preconditioner. 318 \end{memberdesc} 319 320 \begin{memberdesc}[LinearPDE]{DIRECT} 321 the default direct linear solver. 322 \end{memberdesc} 323 324 \begin{memberdesc}[LinearPDE]{CHOLEVSKY} 325 direct solver based on Cholevsky factorization (or similar), see~\Ref{Saad}. The solver will require a symmetric PDE. 326 \end{memberdesc} 327 328 \begin{memberdesc}[LinearPDE]{PCG} 329 preconditioned conjugate gradient method, see~\Ref{WEISS}\index{linear solver!PCG}\index{PCG}. The solver will require a symmetric PDE. 330 \end{memberdesc} 331 332 \begin{memberdesc}[LinearPDE]{GMRES} 333 the GMRES method, see~\Ref{WEISS}\index{linear solver!GMRES}\index{GMRES}. Truncation and restart are controlled by the parameters 334 \var{truncation} and \var{restart} of \method{getSolution}. 335 \end{memberdesc} 336 337 \begin{memberdesc}[LinearPDE]{LUMPING} 338 uses lumping to solve the system of linear equations~\index{linear solver!lumping}\index{lumping}. This solver technique 339 condenses the stiffness matrix to a diagonal matrix so the solution of the linear systems becomes very cheap. It can be used when 340 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 341 but is expect to converge to zero when the mesh gets finer. 342 Lumping does not use the linear system solver library. 343 \end{memberdesc} 344 345 \begin{memberdesc}[LinearPDE]{PRES20} 346 the GMRES method with truncation after five residuals and 347 restart after 20 steps, see~\Ref{WEISS}. 348 \end{memberdesc}[LinearPDE]{CR} 349 350 \begin{memberdesc}[LinearPDE]{CGS} 351 conjugate gradient squared method, see~\Ref{WEISS}. 352 \end{memberdesc} 353 354 \begin{memberdesc}[LinearPDE]{BICGSTAB} 355 stabilized bi-conjugate gradients methods, see~\Ref{WEISS}. 356 \end{memberdesc} 357 358 \begin{memberdesc}[LinearPDE]{SSOR} 359 symmetric successive over-relaxation method, see~\Ref{WEISS}. Typically used as preconditioner but some linear solver libraries support 360 this as a solver. 361 \end{memberdesc} 362 \begin{memberdesc}[LinearPDE]{ILU0} 363 the incomplete LU factorization preconditioner with no fill-in, see~\Ref{Saad}. 364 \end{memberdesc} 365 366 \begin{memberdesc}[LinearPDE]{ILUT} 367 the incomplete LU factorization preconditioner with fill-in, see~\Ref{Saad}. During the LU-factorization element with 368 relative size less then \var{drop_tolerance} are dropped. Moreover, the size of the LU-factorization is restricted to the 369 \var{drop_storage}-fold of the stiffness matrix. \var{drop_tolerance} and \var{drop_storage} are both set in the 370 \method{getSolution} call. 371 \end{memberdesc} 372 373 \begin{memberdesc}[LinearPDE]{JACOBI} 374 the Jacobi preconditioner, see~\Ref{Saad}. 375 \end{memberdesc} 376 377 \begin{memberdesc}[LinearPDE]{AMG} 378 the algebraic--multi grid method, see~\Ref{AMG}. This method can be used as linear solver method but is more robust when used 379 in a preconditioner. 380 \end{memberdesc} 381 382 \begin{memberdesc}[LinearPDE]{RILU} 383 recursive incomplete LU factorization preconditioner, see~\Ref{RILU}. This method is similar to \ILUT but uses smoothing 384 between levels. During the LU-factorization element with 385 relative size less then \var{drop_tolerance} are dropped. Moreover, the size of the LU-factorization is restricted to the 386 \var{drop_storage}-fold of the stiffness matrix. \var{drop_tolerance} and \var{drop_storage} are both set in the 387 \method{getSolution} call. 388 \end{memberdesc} 389 390 \begin{memberdesc}[LinearPDE]{NO_REORDERING} 391 no ordering is used during factorization. 392 \end{memberdesc} 393 394 \begin{memberdesc}[LinearPDE]{MINIMUM_FILL_IN} 395 applies reordering before factorization using a fill-in minimization strategy. You have to check with the particular solver library or 396 linear solver package if this is supported. In any case, it is advisable to apply reordering on the mesh to minimize fill-in. 397 \end{memberdesc} 398 399 \begin{memberdesc}[LinearPDE]{NESTED_DISSECTION} 400 applies reordering before factorization using a nested dissection strategy. You have to check with the particular solver library or 401 linear solver package if this is supported. In any case, it is advisable to apply reordering on the mesh to minimize fill-in. 402 \end{memberdesc} 403 404 \section{The \Poisson Class} 405 The \Poisson class provides an easy way to define and solve the Poisson 406 equation 407 \begin{equation}\label{POISSON.1} 408 -u\hackscore{,ii}=f\; . 409 \end{equation} 410 with homogeneous boundary conditions 411 \begin{equation}\label{POISSON.2} 412 n\hackscore{i}u\hackscore{,i}=0 413 \end{equation} 414 and homogeneous constraints 415 \begin{equation}\label{POISSON.3} 416 u=0 \mbox{ where } q>0 417 \end{equation} 418 $f$ has to be a \Scalar in the \Function and $q$ must be 419 a \Scalar in the \SolutionFS. 420 421 \begin{classdesc}{Poisson}{domain} 422 opens a Poisson equation on the \Domain domain. \Poisson is derived from \LinearPDE. 423 \end{classdesc} 424 \begin{methoddesc}[Poisson]{setValue}{f=escript.Data(),q=escript.Data()} 425 assigns new values to \var{f} and \var{q}. 426 \end{methoddesc} 427 428 \section{The \Helmholtz Class} 429 430 \section{The \Lame Class} 431

## Properties

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