/[escript]/trunk/escript/py_src/linearPDEs.py
ViewVC logotype

Diff of /trunk/escript/py_src/linearPDEs.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 2337 by gross, Thu Mar 26 07:07:42 2009 UTC revision 3402 by gross, Tue Dec 7 07:36:12 2010 UTC
# Line 1  Line 1 
1    # -*- coding: utf-8 -*-
2    
3  ########################################################  ########################################################
4  #  #
5  # Copyright (c) 2003-2008 by University of Queensland  # Copyright (c) 2003-2010 by University of Queensland
6  # Earth Systems Science Computational Center (ESSCC)  # Earth Systems Science Computational Center (ESSCC)
7  # http://www.uq.edu.au/esscc  # http://www.uq.edu.au/esscc
8  #  #
# Line 11  Line 12 
12  #  #
13  ########################################################  ########################################################
14    
15  __copyright__="""Copyright (c) 2003-2008 by University of Queensland  __copyright__="""Copyright (c) 2003-2010 by University of Queensland
16  Earth Systems Science Computational Center (ESSCC)  Earth Systems Science Computational Center (ESSCC)
17  http://www.uq.edu.au/esscc  http://www.uq.edu.au/esscc
18  Primary Business: Queensland, Australia"""  Primary Business: Queensland, Australia"""
19  __license__="""Licensed under the Open Software License version 3.0  __license__="""Licensed under the Open Software License version 3.0
20  http://www.opensource.org/licenses/osl-3.0.php"""  http://www.opensource.org/licenses/osl-3.0.php"""
21  __url__="http://www.uq.edu.au/esscc/escript-finley"  __url__="https://launchpad.net/escript-finley"
22    
23  """  """
24  The module provides an interface to define and solve linear partial  The module provides an interface to define and solve linear partial
25  differential equations (PDEs) and Transport problems within L{escript}.  differential equations (PDEs) and Transport problems within `escript`.
26  L{linearPDEs} does not provide any solver capabilities in itself but hands the  `linearPDEs` does not provide any solver capabilities in itself but hands the
27  PDE over to the PDE solver library defined through the L{Domain<escript.Domain>}  PDE over to the PDE solver library defined through the `Domain`
28  of the PDE. The general interface is provided through the L{LinearPDE} class.  of the PDE. The general interface is provided through the `LinearPDE` class.
29  L{TransportProblem} provides an interface to initial value problems dominated  `TransportProblem` provides an interface to initial value problems dominated
30  by its advective terms.  by its advective terms.
31    
32  @var __author__: name of author  :var __author__: name of author
33  @var __copyright__: copyrights  :var __copyright__: copyrights
34  @var __license__: licence agreement  :var __license__: licence agreement
35  @var __url__: url entry point on documentation  :var __url__: url entry point on documentation
36  @var __version__: version  :var __version__: version
37  @var __date__: date of the version  :var __date__: date of the version
38  """  """
39    
40  import math  import math
41  import escript  import escript
42  import util  import util
43  import numarray  import numpy
44    
45  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
46    
47    
48    class SolverOptions(object):
49        """
50        this class defines the solver options for a linear or non-linear solver.
51        
52        The option also supports the handling of diagnostic informations.
53        
54        Typical usage is
55        
56        ::
57        
58          opts=SolverOptions()
59          print opts
60          opts.resetDiagnostics()
61          u=solver(opts)
62          print "number of iteration steps: =",opts.getDiagnostics("num_iter")
63    
64        :cvar DEFAULT: The default method used to solve the system of linear equations
65        :cvar DIRECT: The direct solver based on LDU factorization
66        :cvar CHOLEVSKY: The direct solver based on LDLt factorization (can only be applied for symmetric PDEs)
67        :cvar PCG: The preconditioned conjugate gradient method (can only be applied for symmetric PDEs)
68        :cvar CR: The conjugate residual method
69        :cvar CGS: The conjugate gradient square method
70        :cvar BICGSTAB: The stabilized Bi-Conjugate Gradient method
71        :cvar TFQMR: Transpose Free Quasi Minimal Residual method
72        :cvar MINRES: Minimum residual method
73        :cvar ILU0: The incomplete LU factorization preconditioner with no fill-in
74        :cvar ILUT: The incomplete LU factorization preconditioner with fill-in
75        :cvar JACOBI: The Jacobi preconditioner
76        :cvar GMRES: The Gram-Schmidt minimum residual method
77        :cvar PRES20: Special GMRES with restart after 20 steps and truncation after 5 residuals
78        :cvar ROWSUM_LUMPING: Matrix lumping using row sum
79        :cvar HRZ_LUMPING: Matrix lumping using the HRZ approach
80        :cvar NO_REORDERING: No matrix reordering allowed
81        :cvar MINIMUM_FILL_IN: Reorder matrix to reduce fill-in during factorization
82        :cvar NESTED_DISSECTION: Reorder matrix to improve load balancing during factorization
83        :cvar PASO: PASO solver package
84        :cvar SCSL: SGI SCSL solver library
85        :cvar MKL: Intel's MKL solver library
86        :cvar UMFPACK: The UMFPACK library
87        :cvar TRILINOS: The TRILINOS parallel solver class library from Sandia National Labs
88        :cvar ITERATIVE: The default iterative solver
89        :cvar AMG: Algebraic Multi Grid
90        :cvar AMLI: Algebraic Multi Level Iteration
91        :cvar REC_ILU: recursive ILU0
92        :cvar RILU: relaxed ILU0
93        :cvar GAUSS_SEIDEL: Gauss-Seidel preconditioner
94        :cvar DEFAULT_REORDERING: the reordering method recommended by the solver
95        :cvar SUPER_LU: the Super_LU solver package
96        :cvar PASTIX: the Pastix direct solver_package
97        :cvar YAIR_SHAPIRA_COARSENING: AMG and AMLI coarsening method by Yair-Shapira
98        :cvar RUGE_STUEBEN_COARSENING: AMG and AMLI coarsening method by Ruge and Stueben
99        :cvar AGGREGATION_COARSENING: AMG and AMLI coarsening using (symmetric) aggregation
100        :cvar STANDARD_COARSENING: AMG and AMLI standard coarsening using mesure of importance of the unknowns
101        :cvar MIN_COARSE_MATRIX_SIZE: minimum size of the coarsest level matrix to use direct solver.
102        :cvar NO_PRECONDITIONER: no preconditioner is applied.
103        """
104        DEFAULT= 0
105        DIRECT= 1
106        CHOLEVSKY= 2
107        PCG= 3
108        CR= 4
109        CGS= 5
110        BICGSTAB= 6
111        ILU0= 8
112        ILUT= 9
113        JACOBI= 10
114        GMRES= 11
115        PRES20= 12
116        LUMPING= 13
117        ROWSUM_LUMPING= 13
118        HRZ_LUMPING= 14
119        NO_REORDERING= 17
120        MINIMUM_FILL_IN= 18
121        NESTED_DISSECTION= 19
122        MKL= 15
123        UMFPACK= 16
124        ITERATIVE= 20
125        PASO= 21
126        AMG= 22
127        REC_ILU = 23
128        TRILINOS = 24
129        NONLINEAR_GMRES = 25
130        TFQMR = 26
131        MINRES = 27
132        GAUSS_SEIDEL=28
133        RILU=29
134        DEFAULT_REORDERING=30
135        SUPER_LU=31
136        PASTIX=32
137        YAIR_SHAPIRA_COARSENING=33
138        RUGE_STUEBEN_COARSENING=34
139        AGGREGATION_COARSENING=35
140        NO_PRECONDITIONER=36
141        MIN_COARSE_MATRIX_SIZE=37
142        AMLI=38
143        STANDARD_COARSENING=39
144    
145        def __init__(self):
146            self.setLevelMax()
147            self.setCoarseningThreshold()
148            self.setSmoother()
149            self.setNumSweeps()
150            self.setNumPreSweeps()
151            self.setNumPostSweeps()
152            self.setTolerance()
153            self.setAbsoluteTolerance()
154            self.setInnerTolerance()
155            self.setDropTolerance()
156            self.setDropStorage()
157            self.setIterMax()
158            self.setInnerIterMax()
159            self.setTruncation()
160            self.setRestart()
161            self.setSymmetry()
162            self.setVerbosity()
163            self.setInnerToleranceAdaption()
164            self.setAcceptanceConvergenceFailure()
165            self.setReordering()
166            self.setPackage()
167            self.setSolverMethod()
168            self.setPreconditioner()
169            self.setCoarsening()
170            self.setMinCoarseMatrixSize()
171            self.setRelaxationFactor()        
172            self.setLocalPreconditionerOff()
173            self.resetDiagnostics(all=True)
174            self.setMinCoarseMatrixSparsity()
175            self.setNumRefinements()
176            self.setNumCoarseMatrixRefinements()
177            self.setUsePanel()
178            self.setUseDirectInterpolation()
179            self.setDiagonalDominanceThreshold()
180            
181    
182        def __str__(self):
183            return self.getSummary()
184        def getSummary(self):
185            """
186            Returns a string reporting the current settings
187            """
188            out="Solver Package: %s"%(self.getName(self.getPackage()))
189            out+="\nVerbosity = %s"%self.isVerbose()
190            out+="\nAccept failed convergence = %s"%self.acceptConvergenceFailure()
191            out+="\nRelative tolerance = %e"%self.getTolerance()
192            out+="\nAbsolute tolerance = %e"%self.getAbsoluteTolerance()
193            out+="\nSymmetric problem = %s"%self.isSymmetric()
194            out+="\nMaximum number of iteration steps = %s"%self.getIterMax()
195            # out+="\nInner tolerance = %e"%self.getInnerTolerance()
196            # out+="\nAdapt innner tolerance = %s"%self.adaptInnerTolerance()
197        
198            if self.getPackage() == self.PASO:
199                out+="\nSolver method = %s"%self.getName(self.getSolverMethod())
200                if self.getSolverMethod() == self.GMRES:
201                    out+="\nTruncation  = %s"%self.getTruncation()
202                    out+="\nRestart  = %s"%self.getRestart()
203                if self.getSolverMethod() == self.AMG:
204                    out+="\nNumber of pre / post sweeps = %s / %s, %s"%(self.getNumPreSweeps(), self.getNumPostSweeps(), self.getNumSweeps())
205                    out+="\nMaximum number of levels = %s"%self.getLevelMax()
206                    out+="\nCoarsening threshold = %e"%self.getCoarseningThreshold()
207                    out+="\nCoarsening method = %s"%self.getName(self.getCoarsening())
208                out+="\nPreconditioner = %s"%self.getName(self.getPreconditioner())
209                out+="\nApply preconditioner locally = %s"%self.useLocalPreconditioner()
210                if self.getPreconditioner() == self.AMG:
211                    out+="\nMaximum number of levels = %s"%self.getLevelMax()
212                    out+="\nCoarsening threshold = %e"%self.getCoarseningThreshold()
213                    out+="\nMinimal sparsity on coarsest level = %e"%(self.getMinCoarseMatrixSparsity(),)
214                    out+="\nSmoother = %s"%self.getName(self.getSmoother())
215                    out+="\nMinimum size of the coarsest level matrix = %e"%self.getCoarseningThreshold()
216                    out+="\nNumber of pre / post sweeps = %s / %s, %s"%(self.getNumPreSweeps(), self.getNumPostSweeps(), self.getNumSweeps())
217                    out+="\nNumber of refinement steps in coarsest level solver = %d"%(self.getNumCoarseMatrixRefinements(),)
218                    out+="\nUse node panel = %s"%(self.usePanel())
219                    out+="\nUse direct interpolation only = %s"%(self.useDirectInterpolation())
220                    out+="\nThreshold for diagonal dominant rows = %s"%(setDiagonalDominanceThreshold(),)
221    
222                if self.getPreconditioner() == self.AMLI:
223                    out+="\nMaximum number of levels = %s"%self.getLevelMax()
224                    out+="\nCoarsening method = %s"%self.getName(self.getCoarsening())
225                    out+="\nCoarsening threshold = %e"%self.getMinCoarseMatrixSize()
226                    out+="\nMinimum size of the coarsest level matrix = %e"%self.getCoarseningThreshold()
227                    out+="\nNumber of pre / post sweeps = %s / %s, %s"%(self.getNumPreSweeps(), self.getNumPostSweeps(), self.getNumSweeps())
228            if self.getPreconditioner() == self.GAUSS_SEIDEL:
229                    out+="\nNumber of sweeps = %s"%self.getNumSweeps()
230                if self.getPreconditioner() == self.ILUT:
231                    out+="\nDrop tolerance = %e"%self.getDropTolerance()
232                    out+="\nStorage increase = %e"%self.getDropStorage()
233                if self.getPreconditioner() == self.RILU:
234                    out+="\nRelaxation factor = %e"%self.getRelaxationFactor()
235            return out
236            
237        def getName(self,key):
238            """
239            returns the name of a given key
240            
241            :param key: a valid key
242            """
243            if key == self.DEFAULT: return "DEFAULT"
244            if key == self.DIRECT: return "DIRECT"
245            if key == self.CHOLEVSKY: return "CHOLEVSKY"
246            if key == self.PCG: return "PCG"
247            if key == self.CR: return "CR"
248            if key == self.CGS: return "CGS"
249            if key == self.BICGSTAB: return "BICGSTAB"
250            if key == self.ILU0: return "ILU0:"
251            if key == self.ILUT: return "ILUT"
252            if key == self.JACOBI: return "JACOBI"
253            if key == self.GMRES: return "GMRES"
254            if key == self.PRES20: return "PRES20"
255            if key == self.ROWSUM_LUMPING: return "ROWSUM_LUMPING"
256            if key == self.HRZ_LUMPING: return "HRZ_LUMPING"
257            if key == self.NO_REORDERING: return "NO_REORDERING"
258            if key == self.MINIMUM_FILL_IN: return "MINIMUM_FILL_IN"
259            if key == self.NESTED_DISSECTION: return "NESTED_DISSECTION"
260            if key == self.MKL: return "MKL"
261            if key == self.UMFPACK: return "UMFPACK"
262            if key == self.ITERATIVE: return "ITERATIVE"
263            if key == self.PASO: return "PASO"
264            if key == self.AMG: return "AMG"
265            if key == self.AMLI: return "AMLI"
266            if key == self.REC_ILU: return "REC_ILU"
267            if key == self.TRILINOS: return "TRILINOS"
268            if key == self.NONLINEAR_GMRES: return "NONLINEAR_GMRES"
269            if key == self.TFQMR: return "TFQMR"
270            if key == self.MINRES: return "MINRES"
271            if key == self.GAUSS_SEIDEL: return "GAUSS_SEIDEL"
272            if key == self.RILU: return "RILU"
273            if key == self.DEFAULT_REORDERING: return "DEFAULT_REORDERING"
274            if key == self.SUPER_LU: return "SUPER_LU"
275            if key == self.PASTIX: return "PASTIX"
276            if key == self.YAIR_SHAPIRA_COARSENING: return "YAIR_SHAPIRA_COARSENING"
277            if key == self.RUGE_STUEBEN_COARSENING: return "RUGE_STUEBEN_COARSENING"
278            if key == self.STANDARD_COARSENING: return "STANDARD_COARSENING"
279            if key == self.AGGREGATION_COARSENING: return "AGGREGATION_COARSENING"
280            if key == self.NO_PRECONDITIONER: return "NO_PRECONDITIONER"
281            
282        def resetDiagnostics(self,all=False):
283            """
284            resets the diagnostics
285            
286            :param all: if ``all`` is ``True`` all diagnostics including accumulative counters are reset.
287            :type all: ``bool``
288            """
289            self.__num_iter=None
290            self.__num_level=None
291            self.__num_inner_iter=None
292            self.__time=None
293            self.__set_up_time=None
294            self.__net_time=None
295            self.__residual_norm=None
296            self.__converged=None
297            self.__preconditioner_size=-1
298            self.__time_step_backtracking_used=None
299            if all:
300                self.__cum_num_inner_iter=0
301                self.__cum_num_iter=0
302                self.__cum_time=0
303                self.__cum_set_up_time=0
304                self.__cum_net_time=0
305    
306        def _updateDiagnostics(self, name, value):
307            """
308            Updates diagnostic information
309            
310            :param name: name of  diagnostic information
311            :type name: ``str`` in the list "num_iter", "num_level", "num_inner_iter", "time", "set_up_time", "net_time", "residual_norm", "converged".
312            :param value: new value of the diagnostic information
313            :note: this function is used by a solver to report diagnostics informations.
314            """
315            if name == "num_iter":
316                self.__num_iter=int(value)
317                self.__cum_num_iter+=self.__num_iter
318            if name == "num_level":
319                self.__num_level=int(value)
320            if name == "num_inner_iter":
321                self.__num_inner_iter=int(value)
322                self.__cum_num_inner_iter+=self.__num_inner_iter
323            if name == "time":
324                self.__time=float(value)
325                self.__cum_time+=self.__time
326            if name == "set_up_time":
327                self.__set_up_time=float(value)
328                self.__cum_set_up_time+=self.__set_up_time
329            if name == "net_time":
330                self.__net_time=float(value)
331                self.__cum_net_time+=self.__net_time
332            if name == "residual_norm":
333                self.__residual_norm=float(value)
334            if name == "converged":
335                self.__converged = (value == True)
336            if name == "time_step_backtracking_used":
337                self.__time_step_backtracking_used = (value == True)
338            if name == "coarse_level_sparsity":
339               self.__coarse_level_sparsity = value
340        if name == "num_coarse_unknowns":
341               self.__num_coarse_unknowns= value
342        def getDiagnostics(self, name):
343            """
344            Returns the diagnostic information ``name``. Possible values are:
345            
346        - "num_iter": the number of iteration steps
347            - "cum_num_iter": the cumulative number of iteration steps
348            - "num_level": the number of level in multi level solver
349            - "num_inner_iter": the number of inner iteration steps
350            - "cum_num_inner_iter": the cumulative number of inner iteration steps
351            - "time": execution time
352            - "cum_time": cumulative execution time
353            - "set_up_time": time to set up of the solver, typically this includes factorization and reordering
354            - "cum_set_up_time": cumulative time to set up of the solver
355            - "net_time": net execution time, excluding setup time for the solver and execution time for preconditioner
356            - "cum_net_time": cumulative net execution time
357            - "preconditioner_size": size of preconditioner [Bytes]
358            - "converged": return True if solution has converged.
359            - "time_step_backtracking_used": returns True if time step back tracking has been used.
360            - "coarse_level_sparsity": returns the sparsity of the matrix on the coarsest level
361            - "num_coarse_unknowns": returns the number of unknowns on the coarsest level
362        
363            
364            :param name: name of diagnostic information to return
365            :type name: ``str`` in the list above.
366            :return: requested value. ``None`` is returned if the value is undefined.
367            :note: If the solver has thrown an exception diagnostic values have an undefined status.
368            """
369            if name == "num_iter": return self.__num_iter
370            elif name == "cum_num_iter": return self.__cum_num_iter
371            elif name == "num_level": return self.__num_level
372            elif name == "num_inner_iter": return self.__num_inner_iter
373            elif name == "cum_num_inner_iter": return self.__cum_num_inner_iter
374            elif name == "time": return self.__time
375            elif name == "cum_time": return self.__cum_time
376            elif name == "set_up_time": return self.__set_up_time
377            elif name == "cum_set_up_time": return self.__cum_set_up_time
378            elif name == "net_time": return self.__net_time
379            elif name == "cum_net_time": return self.__cum_net_time
380            elif name == "residual_norm": return self.__residual_norm
381            elif name == "converged": return self.__converged      
382            elif name == "preconditioner_size": return  self.__preconditioner_size
383            elif name == "time_step_backtracking_used": return  self.__time_step_backtracking_used
384            elif name == "coarse_level_sparsity": return self.__coarse_level_sparsity
385            elif name == "num_coarse_unknowns": return self.__num_coarse_unknowns
386            else:
387                raise ValueError,"unknown diagnostic item %s"%name
388        def hasConverged(self):
389            """
390            Returns ``True`` if the last solver call has been finalized successfully.
391            :note: if an exception has been thrown by the solver the status of this flag is undefined.
392            """
393            return self.getDiagnostics("converged")
394        def setCoarsening(self,method=0):
395            """
396            Sets the key of the coarsening method to be applied in AMG or AMLI
397    
398            :param method: selects the coarsening method .
399            :type method: in {SolverOptions.DEFAULT}, `SolverOptions.YAIR_SHAPIRA_COARSENING`,  `SolverOptions.RUGE_STUEBEN_COARSENING`, `SolverOptions.AGGREGATION_COARSENING`
400            """
401        if method==None: method=0
402            if not method in [self.DEFAULT, self.YAIR_SHAPIRA_COARSENING, self.RUGE_STUEBEN_COARSENING, self.AGGREGATION_COARSENING, self.STANDARD_COARSENING,]:
403                 raise ValueError,"unknown coarsening method %s"%method
404            self.__coarsening=method
405        
406        def getCoarsening(self):
407            """
408            Returns the key of the coarsening algorithm to be applied AMG or AMLI
409    
410            :rtype: in the list `SolverOptions.DEFAULT`, `SolverOptions.YAIR_SHAPIRA_COARSENING`, `SolverOptions.RUGE_STUEBEN_COARSENING`, `SolverOptions.AGGREGATION_COARSENING`
411            """
412            return self.__coarsening
413          
414        def setMinCoarseMatrixSize(self,size=500):
415            """
416            Sets the minumum size of the coarsest level matrix in AMG or AMLI
417    
418            :param size: minumum size of the coarsest level matrix .
419            :type size: positive ``int`` or ``None``
420            """
421        if size==None: size=500
422            size=int(size)
423            if size<0:
424               raise ValueError,"minumum size of the coarsest level matrix must be non-negative."
425            self.__MinCoarseMatrixSize=size
426            
427        def getMinCoarseMatrixSize(self):
428            """
429            Returns the minumum size of the coarsest level matrix in AMG or AMLI
430    
431            :rtype: ``int``
432            """
433            return self.__MinCoarseMatrixSize
434          
435        def setPreconditioner(self, preconditioner=10):
436            """
437            Sets the preconditioner to be used.
438    
439            :param preconditioner: key of the preconditioner to be used.
440            :type preconditioner: in `SolverOptions.ILU0`, `SolverOptions.ILUT`, `SolverOptions.JACOBI`,
441                                        `SolverOptions.AMG`, `SolverOptions.AMLI`, `SolverOptions.REC_ILU`, `SolverOptions.GAUSS_SEIDEL`, `SolverOptions.RILU`,
442                                        `SolverOptions.NO_PRECONDITIONER`
443            :note: Not all packages support all preconditioner. It can be assumed that a package makes a reasonable choice if it encounters an unknown preconditioner.
444            """
445        if preconditioner==None: preconditioner=10
446            if not preconditioner in [  SolverOptions.ILU0, SolverOptions.ILUT, SolverOptions.JACOBI,
447                                        SolverOptions.AMG, SolverOptions.AMLI, SolverOptions.REC_ILU, SolverOptions.GAUSS_SEIDEL, SolverOptions.RILU,
448                                        SolverOptions.NO_PRECONDITIONER] :
449                 raise ValueError,"unknown preconditioner %s"%preconditioner
450            self.__preconditioner=preconditioner    
451        def getPreconditioner(self):
452            """
453            Returns key of the preconditioner to be used.
454    
455            :rtype: in the list `SolverOptions.ILU0`, `SolverOptions.ILUT`, `SolverOptions.JACOBI`, SolverOptions.AMLI,
456                                        `SolverOptions.AMG`, `SolverOptions.REC_ILU`, `SolverOptions.GAUSS_SEIDEL`, `SolverOptions.RILU`,
457                                        `SolverOptions.NO_PRECONDITIONER`
458            """
459            return self.__preconditioner
460        def setSmoother(self, smoother=None):
461            """
462            Sets the smoother to be used.
463    
464            :param smoother: key of the smoother to be used.
465            :type smoother: in `SolverOptions.JACOBI`, `SolverOptions.GAUSS_SEIDEL`
466            :note: Not all packages support all smoothers. It can be assumed that a package makes a reasonable choice if it encounters an unknown smoother.
467            """
468        if smoother==None: smoother=28
469            if not smoother in [ SolverOptions.JACOBI, SolverOptions.GAUSS_SEIDEL ] :
470                 raise ValueError,"unknown smoother %s"%smoother
471            self.__smoother=smoother    
472        def getSmoother(self):
473            """
474            Returns key of the smoother to be used.
475    
476            :rtype: in the list `SolverOptions.JACOBI`, `SolverOptions.GAUSS_SEIDEL`
477            """
478            return self.__smoother  
479        def setSolverMethod(self, method=0):
480            """
481            Sets the solver method to be used. Use ``method``=``SolverOptions.DIRECT`` to indicate that a direct rather than an iterative
482            solver should be used and Use ``method``=``SolverOptions.ITERATIVE`` to indicate that an iterative rather than a direct
483            solver should be used.
484    
485            :param method: key of the solver method to be used.
486            :type method: in `SolverOptions.DEFAULT`, `SolverOptions.DIRECT`, `SolverOptions.CHOLEVSKY`, `SolverOptions.PCG`,
487                            `SolverOptions.CR`, `SolverOptions.CGS`, `SolverOptions.BICGSTAB`,
488                            `SolverOptions.GMRES`, `SolverOptions.PRES20`, `SolverOptions.ROWSUM_LUMPING`, `SolverOptions.HRZ_LUMPING`, `SolverOptions.ITERATIVE`,
489                            `SolverOptions.NONLINEAR_GMRES`, `SolverOptions.TFQMR`, `SolverOptions.MINRES`
490            :note: Not all packages support all solvers. It can be assumed that a package makes a reasonable choice if it encounters an unknown solver method.
491            """
492        if method==None: method=0
493            if not method in [ SolverOptions.DEFAULT, SolverOptions.DIRECT, SolverOptions.CHOLEVSKY, SolverOptions.PCG,
494                               SolverOptions.CR, SolverOptions.CGS, SolverOptions.BICGSTAB,
495                               SolverOptions.GMRES, SolverOptions.PRES20, SolverOptions.ROWSUM_LUMPING, SolverOptions.HRZ_LUMPING,
496                               SolverOptions.ITERATIVE,
497                               SolverOptions.NONLINEAR_GMRES, SolverOptions.TFQMR, SolverOptions.MINRES ]:
498                 raise ValueError,"unknown solver method %s"%method
499            self.__method=method
500        def getSolverMethod(self):
501            """
502            Returns key of the solver method to be used.
503    
504            :rtype: in the list `SolverOptions.DEFAULT`, `SolverOptions.DIRECT`, `SolverOptions.CHOLEVSKY`, `SolverOptions.PCG`,
505                            `SolverOptions.CR`, `SolverOptions.CGS`, `SolverOptions.BICGSTAB`,
506                            `SolverOptions.GMRES`, `SolverOptions.PRES20`, `SolverOptions.ROWSUM_LUMPING`, `SolverOptions.HRZ_LUMPING`, `SolverOptions.ITERATIVE`,
507                            `SolverOptions.NONLINEAR_GMRES`, `SolverOptions.TFQMR`, `SolverOptions.MINRES`
508            """
509            return self.__method
510            
511        def setPackage(self, package=0):
512            """
513            Sets the solver package to be used as a solver.  
514    
515            :param package: key of the solver package to be used.
516            :type package: in `SolverOptions.DEFAULT`, `SolverOptions.PASO`, `SolverOptions.SUPER_LU`, `SolverOptions.PASTIX`, `SolverOptions.MKL`, `SolverOptions.UMFPACK`, `SolverOptions.TRILINOS`
517            :note: Not all packages are support on all implementation. An exception may be thrown on some platforms if a particular is requested.
518            """
519        if package==None: package=0
520            if not package in [SolverOptions.DEFAULT, SolverOptions.PASO, SolverOptions.SUPER_LU, SolverOptions.PASTIX, SolverOptions.MKL, SolverOptions.UMFPACK, SolverOptions.TRILINOS]:
521                 raise ValueError,"unknown solver package %s"%package
522            self.__package=package
523        def getPackage(self):
524            """
525            Returns the solver package key
526    
527            :rtype: in the list `SolverOptions.DEFAULT`, `SolverOptions.PASO`, `SolverOptions.SUPER_LU`, `SolverOptions.PASTIX`, `SolverOptions.MKL`, `SolverOptions.UMFPACK`, `SolverOptions.TRILINOS`
528            """
529            return self.__package
530        def setReordering(self,ordering=30):
531            """
532            Sets the key of the reordering method to be applied if supported by the solver. Some direct solvers support reordering
533            to optimize compute time and storage use during elimination.
534    
535            :param ordering: selects the reordering strategy.
536            :type ordering: in 'SolverOptions.NO_REORDERING', 'SolverOptions.MINIMUM_FILL_IN', 'SolverOptions.NESTED_DISSECTION', 'SolverOptions.DEFAULT_REORDERING'
537            """
538            if not ordering in [self.NO_REORDERING, self.MINIMUM_FILL_IN, self.NESTED_DISSECTION, self.DEFAULT_REORDERING]:
539                 raise ValueError,"unknown reordering strategy %s"%ordering
540            self.__reordering=ordering
541        def getReordering(self):
542            """
543            Returns the key of the reordering method to be applied if supported by the solver.
544    
545            :rtype: ordering in 'SolverOptions.NO_REORDERING', 'SolverOptions.MINIMUM_FILL_IN', 'SolverOptions.NESTED_DISSECTION', 'SolverOptions.DEFAULT_REORDERING'
546            """
547            return self.__reordering
548        def setRestart(self,restart=None):
549            """
550            Sets the number of iterations steps after which GMRES is performing a restart.
551    
552            :param restart: number of iteration steps after which to perform a restart. If equal to ``None`` no restart is performed.
553            :type restart: ``int`` or ``None``
554            """
555            if restart == None:
556                self.__restart=restart
557            else:
558                restart=int(restart)
559                if restart<1:
560                    raise ValueError,"restart must be positive."
561                self.__restart=restart
562            
563        def getRestart(self):
564            """
565            Returns the number of iterations steps after which GMRES is performing a restart.
566            If ``None`` is returned no restart is performed.
567    
568            :rtype: ``int`` or ``None``
569            """
570            if self.__restart < 0:
571                return None
572            else:
573                return self.__restart
574        def _getRestartForC(self):
575            r=self.getRestart()
576            if r == None:
577                return -1
578                else:
579                return r
580      
581        def setDiagonalDominanceThreshold(self,value=0.5):
582            """
583            Sets the threshold for diagonal dominant rows which are eliminated during AMG coarsening.
584    
585         :param value: threshold
586         :type value: ``float``
587            """
588            value=float(value)
589            if value<0 or value>1.:
590           raise ValueError,"Diagonal dominance threshold must be between 0 and 1."
591        self.__diagonal_dominance_threshold=value
592            
593        def getDiagonalDominanceThreshold(self):
594            """
595            Returns the threshold for diagonal dominant rows which are eliminated during AMG coarsening.
596    
597            :rtype: ``float``
598            """
599            return self.__diagonal_dominance_threshold
600            
601        def setTruncation(self,truncation=20):
602            """
603            Sets the number of residuals in GMRES to be stored for orthogonalization.  The more residuals are stored
604            the faster GMRES converged but
605    
606            :param truncation: truncation
607            :type truncation: ``int``
608            """
609            truncation=int(truncation)
610            if truncation<1:
611               raise ValueError,"truncation must be positive."
612            self.__truncation=truncation
613        def getTruncation(self):
614            """
615            Returns the number of residuals in GMRES to be stored for orthogonalization
616    
617            :rtype: ``int``
618            """
619            return self.__truncation
620        def setInnerIterMax(self,iter_max=10):
621            """
622            Sets the maximum number of iteration steps for the inner iteration.
623    
624            :param iter_max: maximum number of inner iterations
625            :type iter_max: ``int``
626            """
627            iter_max=int(iter_max)
628            if iter_max<1:
629               raise ValueError,"maximum number of inner iteration must be positive."
630            self.__inner_iter_max=iter_max
631        def getInnerIterMax(self):
632            """
633            Returns maximum number of inner iteration steps
634    
635            :rtype: ``int``
636            """
637            return self.__inner_iter_max
638        def setIterMax(self,iter_max=100000):
639            """
640            Sets the maximum number of iteration steps
641    
642            :param iter_max: maximum number of iteration steps
643            :type iter_max: ``int``
644            """
645            iter_max=int(iter_max)
646            if iter_max<1:
647               raise ValueError,"maximum number of iteration steps must be positive."
648            self.__iter_max=iter_max
649        def getIterMax(self):
650            """
651            Returns maximum number of iteration steps
652    
653            :rtype: ``int``
654            """
655            return self.__iter_max
656        def setLevelMax(self,level_max=100):
657            """
658            Sets the maximum number of coarsening levels to be used in an algebraic multi level solver or preconditioner
659    
660            :param level_max: maximum number of levels
661            :type level_max: ``int``
662            """
663            level_max=int(level_max)
664            if level_max<0:
665               raise ValueError,"maximum number of coarsening levels must be non-negative."
666            self.__level_max=level_max
667        def getLevelMax(self):
668            """
669            Returns the maximum number of coarsening levels to be used in an algebraic multi level solver or preconditioner
670    
671            :rtype: ``int``
672            """
673            return self.__level_max
674        def setCoarseningThreshold(self,theta=0.25):
675            """
676            Sets the threshold for coarsening in the algebraic multi level solver or preconditioner
677    
678            :param theta: threshold for coarsening
679            :type theta: positive ``float``
680            """
681            theta=float(theta)
682            if theta<0 or theta>1:
683               raise ValueError,"threshold must be non-negative and less or equal 1."
684            self.__coarsening_threshold=theta
685        def getCoarseningThreshold(self):
686            """
687            Returns the threshold for coarsening in the algebraic multi level solver or preconditioner
688    
689            :rtype: ``float``
690            """
691            return self.__coarsening_threshold
692        def setNumSweeps(self,sweeps=1):
693            """
694            Sets the number of sweeps in a Jacobi or Gauss-Seidel/SOR preconditioner.
695    
696            :param sweeps: number of sweeps
697            :type sweeps: positive ``int``
698            """
699            sweeps=int(sweeps)
700            if sweeps<1:
701               raise ValueError,"number of sweeps must be positive."
702            self.__sweeps=sweeps
703        def getNumSweeps(self):
704            """
705            Returns the number of sweeps in a Jacobi or Gauss-Seidel/SOR preconditioner.
706    
707            :rtype: ``int``
708            """
709            return self.__sweeps
710        def setNumPreSweeps(self,sweeps=2):
711            """
712            Sets the number of sweeps in the pre-smoothing step of a multi level solver or preconditioner
713    
714            :param sweeps: number of sweeps
715            :type sweeps: positive ``int``
716            """
717            sweeps=int(sweeps)
718            if sweeps<1:
719               raise ValueError,"number of sweeps must be positive."
720            self.__pre_sweeps=sweeps
721        def getNumPreSweeps(self):
722            """
723            Returns he number of sweeps in the pre-smoothing step of a multi level solver or preconditioner
724    
725            :rtype: ``int``
726            """
727            return self.__pre_sweeps
728        def setNumPostSweeps(self,sweeps=2):
729            """
730            Sets the number of sweeps in the post-smoothing step of a multi level solver or preconditioner
731    
732            :param sweeps: number of sweeps
733            :type sweeps: positive ``int``
734            """
735            sweeps=int(sweeps)
736            if sweeps<1:
737               raise ValueError,"number of sweeps must be positive."
738            self.__post_sweeps=sweeps
739        def getNumPostSweeps(self):
740            """
741            Returns he number of sweeps in the post-smoothing step of a multi level solver or preconditioner
742    
743            :rtype: ``int``
744            """
745            return self.__post_sweeps
746    
747        def setTolerance(self,rtol=1.e-8):
748            """
749            Sets the relative tolerance for the solver
750    
751            :param rtol: relative tolerance
752            :type rtol: non-negative ``float``
753            """
754            rtol=float(rtol)
755            if rtol<0 or rtol>1:
756               raise ValueError,"tolerance must be non-negative and less or equal 1."
757            self.__tolerance=rtol
758        def getTolerance(self):
759            """
760            Returns the relative tolerance for the solver
761    
762            :rtype: ``float``
763            """
764            return self.__tolerance
765        def setAbsoluteTolerance(self,atol=0.):
766            """
767            Sets the absolute tolerance for the solver
768    
769            :param atol:  absolute tolerance
770            :type atol: non-negative ``float``
771            """
772            atol=float(atol)
773            if atol<0:
774               raise ValueError,"tolerance must be non-negative."
775            self.__absolute_tolerance=atol
776        def getAbsoluteTolerance(self):
777            """
778            Returns the absolute tolerance for the solver
779    
780            :rtype: ``float``
781            """
782            return self.__absolute_tolerance
783    
784        def setInnerTolerance(self,rtol=0.9):
785            """
786             Sets the relative tolerance for an inner iteration scheme for instance on the coarsest level in a multi-level scheme.
787    
788            :param rtol: inner relative tolerance
789            :type rtol: positive ``float``
790            """
791            rtol=float(rtol)
792            if rtol<=0 or rtol>1:
793                raise ValueError,"tolerance must be positive and less or equal 1."
794            self.__inner_tolerance=rtol
795        def getInnerTolerance(self):
796            """
797            Returns the relative tolerance for an inner iteration scheme
798    
799            :rtype: ``float``
800            """
801            return self.__inner_tolerance
802        def setDropTolerance(self,drop_tol=0.01):
803            """
804            Sets the relative drop tolerance in ILUT
805    
806            :param drop_tol: drop tolerance
807            :type drop_tol: positive ``float``
808            """
809            drop_tol=float(drop_tol)
810            if drop_tol<=0 or drop_tol>1:
811                raise ValueError,"drop tolerance must be positive and less or equal 1."
812            self.__drop_tolerance=drop_tol
813        def getDropTolerance(self):
814            """
815            Returns the relative drop tolerance in ILUT
816    
817            :rtype: ``float``
818            """
819            return self.__drop_tolerance
820        def setDropStorage(self,storage=2.):
821            """
822            Sets the maximum allowed increase in storage for ILUT. ``storage`` =2 would mean that
823            a doubling of the storage needed for the coefficient matrix is allowed in the ILUT factorization.
824    
825            :param storage: allowed storage increase
826            :type storage: ``float``
827            """
828            storage=float(storage)
829            if storage<1:
830                raise ValueError,"allowed storage increase must be greater or equal to 1."
831            self.__drop_storage=storage
832        def getDropStorage(self):
833        
834            """
835            Returns the maximum allowed increase in storage for ILUT
836    
837            :rtype: ``float``
838            """
839            return self.__drop_storage
840        def setRelaxationFactor(self,factor=0.3):
841            """
842            Sets the relaxation factor used to add dropped elements in RILU to the main diagonal.
843    
844            :param factor: relaxation factor
845            :type factor: ``float``
846            :note: RILU with a relaxation factor 0 is identical to ILU0
847            """
848            factor=float(factor)
849            if factor<0:
850                raise ValueError,"relaxation factor must be non-negative."
851            self.__relaxation=factor
852        def getRelaxationFactor(self):
853        
854            """
855            Returns the relaxation factor used to add dropped elements in RILU to the main diagonal.
856    
857            :rtype: ``float``
858            """
859            return self.__relaxation
860        def isSymmetric(self):
861            """
862            Checks if symmetry of the  coefficient matrix is indicated.
863    
864            :return: True if a symmetric PDE is indicated, False otherwise
865            :rtype: ``bool``
866            """
867            return self.__symmetric
868        def setSymmetryOn(self):
869            """
870            Sets the symmetry flag to indicate that the coefficient matrix is symmetric.
871            """
872            self.__symmetric=True
873        def setSymmetryOff(self):
874            """
875            Clears the symmetry flag for the coefficient matrix.
876            """
877            self.__symmetric=False
878        def setSymmetry(self,flag=False):
879            """
880            Sets the symmetry flag for the coefficient matrix to ``flag``.
881    
882            :param flag: If True, the symmetry flag is set otherwise reset.
883            :type flag: ``bool``
884            """
885            if flag:
886                self.setSymmetryOn()
887            else:
888                self.setSymmetryOff()
889        def isVerbose(self):
890            """
891            Returns ``True`` if the solver is expected to be verbose.
892    
893            :return: True if verbosity of switched on.
894            :rtype: ``bool``
895            """
896            return self.__verbose
897    
898        def setVerbosityOn(self):
899            """
900            Switches the verbosity of the solver on.
901            """
902            self.__verbose=True
903        def setVerbosityOff(self):
904            """
905            Switches the verbosity of the solver off.
906            """
907            self.__verbose=False
908        def setVerbosity(self,verbose=False):
909            """
910            Sets the verbosity flag for the solver to ``flag``.
911    
912            :param verbose: If ``True``, the verbosity of the solver is switched on.
913            :type verbose: ``bool``
914            """
915            if verbose:
916                self.setVerbosityOn()
917            else:
918                self.setVerbosityOff()
919            
920        def adaptInnerTolerance(self):
921            """
922            Returns ``True`` if the tolerance of the inner solver is selected automatically.
923            Otherwise the inner tolerance set by `setInnerTolerance` is used.
924    
925            :return: ``True`` if inner tolerance adaption is chosen.
926            :rtype: ``bool``
927            """
928            return self.__adapt_inner_tolerance
929    
930        def setInnerToleranceAdaptionOn(self):
931            """
932            Switches the automatic selection of inner tolerance on
933            """
934            self.__adapt_inner_tolerance=True
935        def setInnerToleranceAdaptionOff(self):
936            """
937            Switches the automatic selection of inner tolerance off.
938            """
939            self.__adapt_inner_tolerance=False
940        def setInnerToleranceAdaption(self,adapt=True):
941            """
942            Sets the flag to indicate automatic selection of the inner tolerance.
943    
944            :param adapt: If ``True``, the inner tolerance is selected automatically.
945            :type adapt: ``bool``
946            """
947            if adapt:
948                self.setInnerToleranceAdaptionOn()
949            else:
950                self.setInnerToleranceAdaptionOff()
951    
952        def acceptConvergenceFailure(self):
953            """
954            Returns ``True`` if a failure to meet the stopping criteria within the
955            given number of iteration steps is not raising in exception. This is useful
956            if a solver is used in a non-linear context where the non-linear solver can
957            continue even if the returned the solution does not necessarily meet the
958            stopping criteria. One can use the `hasConverged` method to check if the
959            last call to the solver was successful.
960    
961            :return: ``True`` if a failure to achieve convergence is accepted.
962            :rtype: ``bool``
963            """
964            return self.__accept_convergence_failure
965    
966        def setAcceptanceConvergenceFailureOn(self):
967            """
968            Switches the acceptance of a failure of convergence on
969            """
970            self.__accept_convergence_failure=True
971        def setAcceptanceConvergenceFailureOff(self):
972            """
973            Switches the acceptance of a failure of convergence off.
974            """
975            self.__accept_convergence_failure=False
976        def setAcceptanceConvergenceFailure(self,accept=False):
977            """
978            Sets the flag to indicate the acceptance of a failure of convergence.
979    
980            :param accept: If ``True``, any failure to achieve convergence is accepted.
981            :type accept: ``bool``
982            """
983            if accept:
984                self.setAcceptanceConvergenceFailureOn()
985            else:
986                self.setAcceptanceConvergenceFailureOff()
987    
988        def useLocalPreconditioner(self):
989            """
990            Returns ``True`` if the preconditoner is applied locally on each MPI. This reducess communication costs
991            and speeds up the application of the preconditioner but at the costs of more iteration steps. This can be an
992            advantage on clusters with slower interconnects.
993    
994            :return: ``True`` if local preconditioning is applied
995            :rtype: ``bool``
996            """
997            return self.__use_local_preconditioner
998    
999        def setLocalPreconditionerOn(self):
1000            """
1001            Sets the flag to use  local preconditioning to on
1002            """
1003            self.__use_local_preconditioner=True
1004        def setLocalPreconditionerOff(self):
1005            """
1006            Sets the flag to use  local preconditioning to off
1007            """
1008            self.__use_local_preconditioner=False
1009            
1010        def setLocalPreconditioner(self,use=False):
1011            """
1012            Sets the flag to use  local preconditioning
1013    
1014            :param use: If ``True``, local proconditioning on each MPI rank is applied
1015            :type use: ``bool``
1016            """
1017            if use:
1018                self.setUseLocalPreconditionerOn()
1019            else:
1020                self.setUseLocalPreconditionerOff()
1021                
1022        def setMinCoarseMatrixSparsity(self,sparsity=0.05):
1023           """
1024           Sets the minimum sparsity on the coarsest level. Typically
1025           a direct solver is used when the sparsity becomes bigger than
1026           the set limit.
1027          
1028           :param sparsity: minimual sparsity
1029           :type sparsity: ``float``
1030           """
1031           sparsity=float(sparsity)
1032           if factor<0:
1033         raise ValueError,"sparsity must be non-negative."
1034           if factor>1:
1035               raise ValueError,"sparsity must be less than 1."
1036           self.__min_sparsity=factor
1037        def setMinCoarseMatrixSparsity(self,sparsity=0.05):
1038           """
1039           Sets the minimum sparsity on the coarsest level. Typically
1040           a direct solver is used when the sparsity becomes bigger than
1041           the set limit.
1042          
1043           :param sparsity: minimual sparsity
1044           :type sparsity: ``float``
1045           """
1046           sparsity=float(sparsity)
1047           if sparsity<0:
1048         raise ValueError,"sparsity must be non-negative."
1049           if sparsity>1:
1050              raise ValueError,"sparsity must be less than 1."
1051           self.__min_sparsity=sparsity
1052    
1053        def getMinCoarseMatrixSparsity(self):
1054          """
1055          Returns the minimum sparsity on the coarsest level. Typically
1056          a direct solver is used when the sparsity becomes bigger than
1057          the set limit.
1058      
1059          :return: minimual sparsity
1060          :rtype: ``float``
1061          """
1062          return self.__min_sparsity
1063    
1064        def setNumRefinements(self,refinements=2):
1065          """
1066          Sets the number of refinement steps to refine the solution when a direct solver is applied.
1067      
1068          :param refinements: number of refinements
1069          :type refinements: non-negative ``int``
1070          """
1071          refinements=int(refinements)
1072          if refinements<0:
1073         raise ValueError,"number of refinements must be non-negative."
1074          self.__refinements=refinements
1075      
1076        def getNumRefinements(self):
1077           """
1078           Returns the number of refinement steps to refine the solution when a direct solver is applied.
1079          
1080           :rtype: non-negative ``int``
1081           """
1082           return self.__refinements
1083    
1084        def setNumCoarseMatrixRefinements(self,refinements=2):
1085          """
1086          Sets the number of refinement steps to refine the solution on the coarset level when a direct solver is applied.
1087      
1088          :param refinements: number of refinements
1089          :type refinements: non-negative ``int``
1090          """
1091          refinements=int(refinements)
1092          if refinements<0:
1093         raise ValueError,"number of refinements must be non-negative."
1094          self.__coarse_refinements=refinements
1095      
1096        def getNumCoarseMatrixRefinements(self):
1097          """
1098          Returns the number of resfinement steps to refine the solution on the coarset level when a direct solver is applied.
1099          
1100          :rtype: non-negative ``int``
1101          """
1102          return self.__coarse_refinements
1103    
1104        def usePanel(self):
1105            """
1106            Returns ``True`` if a panel is used to serach for unknown in the AMG coarsening, The panel approach is normally faster
1107            but can lead to larger coarse level systems.
1108            
1109            :return: ``True`` if a panel is used to find unknowns in AMG coarsening
1110            :rtype: ``bool``
1111            """
1112            return self.__use_panel
1113    
1114        def setUsePanelOn(self):
1115            """
1116            Sets the flag to use a panel to find unknowns in AMG coarsening
1117            """
1118            self.__use_panel=True
1119            
1120        def setUsePanelOff(self):
1121            """
1122            Sets the flag to use a panel to find unknowns in AMG coarsening to off
1123            """
1124            self.__use_panel=False
1125            
1126        def setUsePanel(self,use=True):
1127            """
1128            Sets the flag to use  a panel to find unknowns in AMG coarsening
1129    
1130            :param use: If ``True``,a panel is used to find unknowns in AMG coarsening
1131            :type use: ``bool``
1132            """
1133            if use:
1134                self.setUsePanelOn()
1135            else:
1136                self.setUsePanelOff()
1137                
1138        def useDirectInterpolation(self):
1139            """
1140            Returns ``True`` if direct interpolation is used in AMG.
1141    
1142        :return: ``True`` if direct interpolation is used in AMG.
1143            :rtype: ``bool``
1144            """
1145            return self.__use_direct_interpolation
1146    
1147        def setUseDirectInterpolationOn(self):
1148            """
1149            Sets the flag to use direct interpolation in AMG
1150            """
1151            self.__use_direct_interpolation=True
1152            
1153        def setUseDirectInterpolationOff(self):
1154            """
1155            Sets the flag to use direct interpolation in AMG
1156            """
1157            self.__use_direct_interpolation=False
1158            
1159        def setUseDirectInterpolation(self,use=False):
1160            """
1161            Sets the flag to use direct interpolation in AMG
1162    
1163        :param use: If ``True``, direct interpolation in AMG
1164        :type use: ``bool``
1165            """
1166            if use:
1167                self.setUseDirectInterpolationOn()
1168            else:
1169                self.setUseDirectInterpolationOff()
1170    
1171    
1172  class IllegalCoefficient(ValueError):  class IllegalCoefficient(ValueError):
1173     """     """
1174     Exception that is raised if an illegal coefficient of the general or     Exception that is raised if an illegal coefficient of the general or
# Line 73  class PDECoef(object): Line 1198  class PDECoef(object):
1198      """      """
1199      A class for describing a PDE coefficient.      A class for describing a PDE coefficient.
1200    
1201      @cvar INTERIOR: indicator that coefficient is defined on the interior of      :cvar INTERIOR: indicator that coefficient is defined on the interior of
1202                      the PDE domain                      the PDE domain
1203      @cvar BOUNDARY: indicator that coefficient is defined on the boundary of      :cvar BOUNDARY: indicator that coefficient is defined on the boundary of
1204                      the PDE domain                      the PDE domain
1205      @cvar CONTACT: indicator that coefficient is defined on the contact region      :cvar CONTACT: indicator that coefficient is defined on the contact region
1206                     within the PDE domain                     within the PDE domain
1207      @cvar INTERIOR_REDUCED: indicator that coefficient is defined on the      :cvar INTERIOR_REDUCED: indicator that coefficient is defined on the
1208                              interior of the PDE domain using a reduced                              interior of the PDE domain using a reduced
1209                              integration order                              integration order
1210      @cvar BOUNDARY_REDUCED: indicator that coefficient is defined on the      :cvar BOUNDARY_REDUCED: indicator that coefficient is defined on the
1211                              boundary of the PDE domain using a reduced                              boundary of the PDE domain using a reduced
1212                              integration order                              integration order
1213      @cvar CONTACT_REDUCED: indicator that coefficient is defined on the contact      :cvar CONTACT_REDUCED: indicator that coefficient is defined on the contact
1214                             region within the PDE domain using a reduced                             region within the PDE domain using a reduced
1215                             integration order                             integration order
1216      @cvar SOLUTION: indicator that coefficient is defined through a solution of      :cvar SOLUTION: indicator that coefficient is defined through a solution of
1217                      the PDE                      the PDE
1218      @cvar REDUCED: indicator that coefficient is defined through a reduced      :cvar REDUCED: indicator that coefficient is defined through a reduced
1219                     solution of the PDE                     solution of the PDE
1220      @cvar BY_EQUATION: indicator that the dimension of the coefficient shape is      :cvar BY_EQUATION: indicator that the dimension of the coefficient shape is
1221                         defined by the number of PDE equations                         defined by the number of PDE equations
1222      @cvar BY_SOLUTION: indicator that the dimension of the coefficient shape is      :cvar BY_SOLUTION: indicator that the dimension of the coefficient shape is
1223                         defined by the number of PDE solutions                         defined by the number of PDE solutions
1224      @cvar BY_DIM: indicator that the dimension of the coefficient shape is      :cvar BY_DIM: indicator that the dimension of the coefficient shape is
1225                    defined by the spatial dimension                    defined by the spatial dimension
1226      @cvar OPERATOR: indicator that the the coefficient alters the operator of      :cvar OPERATOR: indicator that the the coefficient alters the operator of
1227                      the PDE                      the PDE
1228      @cvar RIGHTHANDSIDE: indicator that the the coefficient alters the right      :cvar RIGHTHANDSIDE: indicator that the the coefficient alters the right
1229                           hand side of the PDE                           hand side of the PDE
1230      @cvar BOTH: indicator that the the coefficient alters the operator as well      :cvar BOTH: indicator that the the coefficient alters the operator as well
1231                  as the right hand side of the PDE                  as the right hand side of the PDE
1232    
1233      """      """
# Line 125  class PDECoef(object): Line 1250  class PDECoef(object):
1250         """         """
1251         Initialises a PDE coefficient type.         Initialises a PDE coefficient type.
1252    
1253         @param where: describes where the coefficient lives         :param where: describes where the coefficient lives
1254         @type where: one of L{INTERIOR}, L{BOUNDARY}, L{CONTACT}, L{SOLUTION},         :type where: one of `INTERIOR`, `BOUNDARY`, `CONTACT`, `SOLUTION`,
1255                      L{REDUCED}, L{INTERIOR_REDUCED}, L{BOUNDARY_REDUCED},                      `REDUCED`, `INTERIOR_REDUCED`, `BOUNDARY_REDUCED`,
1256                      L{CONTACT_REDUCED}                      `CONTACT_REDUCED`
1257         @param pattern: describes the shape of the coefficient and how the shape         :param pattern: describes the shape of the coefficient and how the shape
1258                         is built for a given spatial dimension and numbers of                         is built for a given spatial dimension and numbers of
1259                         equations and solutions in then PDE. For instance,                         equations and solutions in then PDE. For instance,
1260                         (L{BY_EQUATION},L{BY_SOLUTION},L{BY_DIM}) describes a                         (`BY_EQUATION`,`BY_SOLUTION`,`BY_DIM`) describes a
1261                         rank 3 coefficient which is instantiated as shape (3,2,2)                         rank 3 coefficient which is instantiated as shape (3,2,2)
1262                         in case of three equations and two solution components                         in case of three equations and two solution components
1263                         on a 2-dimensional domain. In the case of single equation                         on a 2-dimensional domain. In the case of single equation
1264                         and a single solution component the shape components                         and a single solution component the shape components
1265                         marked by L{BY_EQUATION} or L{BY_SOLUTION} are dropped.                         marked by `BY_EQUATION` or `BY_SOLUTION` are dropped.
1266                         In this case the example would be read as (2,).                         In this case the example would be read as (2,).
1267         @type pattern: C{tuple} of L{BY_EQUATION}, L{BY_SOLUTION}, L{BY_DIM}         :type pattern: ``tuple`` of `BY_EQUATION`, `BY_SOLUTION`, `BY_DIM`
1268         @param altering: indicates what part of the PDE is altered if the         :param altering: indicates what part of the PDE is altered if the
1269                          coefficient is altered                          coefficient is altered
1270         @type altering: one of L{OPERATOR}, L{RIGHTHANDSIDE}, L{BOTH}         :type altering: one of `OPERATOR`, `RIGHTHANDSIDE`, `BOTH`
1271         """         """
1272         super(PDECoef, self).__init__()         super(PDECoef, self).__init__()
1273         self.what=where         self.what=where
# Line 158  class PDECoef(object): Line 1283  class PDECoef(object):
1283    
1284      def getFunctionSpace(self,domain,reducedEquationOrder=False,reducedSolutionOrder=False):      def getFunctionSpace(self,domain,reducedEquationOrder=False,reducedSolutionOrder=False):
1285         """         """
1286         Returns the L{FunctionSpace<escript.FunctionSpace>} of the coefficient.         Returns the `FunctionSpace` of the coefficient.
1287    
1288         @param domain: domain on which the PDE uses the coefficient         :param domain: domain on which the PDE uses the coefficient
1289         @type domain: L{Domain<escript.Domain>}         :type domain: `Domain`
1290         @param reducedEquationOrder: True to indicate that reduced order is used         :param reducedEquationOrder: True to indicate that reduced order is used
1291                                      to represent the equation                                      to represent the equation
1292         @type reducedEquationOrder: C{bool}         :type reducedEquationOrder: ``bool``
1293         @param reducedSolutionOrder: True to indicate that reduced order is used         :param reducedSolutionOrder: True to indicate that reduced order is used
1294                                      to represent the solution                                      to represent the solution
1295         @type reducedSolutionOrder: C{bool}         :type reducedSolutionOrder: ``bool``
1296         @return: L{FunctionSpace<escript.FunctionSpace>} of the coefficient         :return: `FunctionSpace` of the coefficient
1297         @rtype: L{FunctionSpace<escript.FunctionSpace>}         :rtype: `FunctionSpace`
1298         """         """
1299         if self.what==self.INTERIOR:         if self.what==self.INTERIOR:
1300              return escript.Function(domain)              return escript.Function(domain)
# Line 195  class PDECoef(object): Line 1320  class PDECoef(object):
1320         """         """
1321         Returns the value of the coefficient.         Returns the value of the coefficient.
1322    
1323         @return: value of the coefficient         :return: value of the coefficient
1324         @rtype: L{Data<escript.Data>}         :rtype: `Data`
1325         """         """
1326         return self.value         return self.value
1327    
# Line 204  class PDECoef(object): Line 1329  class PDECoef(object):
1329         """         """
1330         Sets the value of the coefficient to a new value.         Sets the value of the coefficient to a new value.
1331    
1332         @param domain: domain on which the PDE uses the coefficient         :param domain: domain on which the PDE uses the coefficient
1333         @type domain: L{Domain<escript.Domain>}         :type domain: `Domain`
1334         @param numEquations: number of equations of the PDE         :param numEquations: number of equations of the PDE
1335         @type numEquations: C{int}         :type numEquations: ``int``
1336         @param numSolutions: number of components of the PDE solution         :param numSolutions: number of components of the PDE solution
1337         @type numSolutions: C{int}         :type numSolutions: ``int``
1338         @param reducedEquationOrder: True to indicate that reduced order is used         :param reducedEquationOrder: True to indicate that reduced order is used
1339                                      to represent the equation                                      to represent the equation
1340         @type reducedEquationOrder: C{bool}         :type reducedEquationOrder: ``bool``
1341         @param reducedSolutionOrder: True to indicate that reduced order is used         :param reducedSolutionOrder: True to indicate that reduced order is used
1342                                      to represent the solution                                      to represent the solution
1343         @type reducedSolutionOrder: C{bool}         :type reducedSolutionOrder: ``bool``
1344         @param newValue: number of components of the PDE solution         :param newValue: number of components of the PDE solution
1345         @type newValue: any object that can be converted into a         :type newValue: any object that can be converted into a
1346                         L{Data<escript.Data>} object with the appropriate shape                         `Data` object with the appropriate shape
1347                         and L{FunctionSpace<escript.FunctionSpace>}                         and `FunctionSpace`
1348         @raise IllegalCoefficientValue: if the shape of the assigned value does         :raise IllegalCoefficientValue: if the shape of the assigned value does
1349                                         not match the shape of the coefficient                                         not match the shape of the coefficient
1350         @raise IllegalCoefficientFunctionSpace: if unable to interpolate value         :raise IllegalCoefficientFunctionSpace: if unable to interpolate value
1351                                                 to appropriate function space                                                 to appropriate function space
1352         """         """
1353         if newValue==None:         if newValue==None:
# Line 245  class PDECoef(object): Line 1370  class PDECoef(object):
1370          """          """
1371          Checks if the coefficient alters the operator of the PDE.          Checks if the coefficient alters the operator of the PDE.
1372    
1373          @return: True if the operator of the PDE is changed when the          :return: True if the operator of the PDE is changed when the
1374                   coefficient is changed                   coefficient is changed
1375          @rtype: C{bool}          :rtype: ``bool``
1376          """          """
1377          if self.altering==self.OPERATOR or self.altering==self.BOTH:          if self.altering==self.OPERATOR or self.altering==self.BOTH:
1378              return not None              return not None
# Line 258  class PDECoef(object): Line 1383  class PDECoef(object):
1383          """          """
1384          Checks if the coefficient alters the right hand side of the PDE.          Checks if the coefficient alters the right hand side of the PDE.
1385    
1386          @rtype: C{bool}          :rtype: ``bool``
1387          @return: True if the right hand side of the PDE is changed when the          :return: True if the right hand side of the PDE is changed when the
1388                   coefficient is changed, C{None} otherwise.                   coefficient is changed, ``None`` otherwise.
1389          """          """
1390          if self.altering==self.RIGHTHANDSIDE or self.altering==self.BOTH:          if self.altering==self.RIGHTHANDSIDE or self.altering==self.BOTH:
1391              return not None              return not None
# Line 272  class PDECoef(object): Line 1397  class PDECoef(object):
1397         Tries to estimate the number of equations and number of solutions if         Tries to estimate the number of equations and number of solutions if
1398         the coefficient has the given shape.         the coefficient has the given shape.
1399    
1400         @param domain: domain on which the PDE uses the coefficient         :param domain: domain on which the PDE uses the coefficient
1401         @type domain: L{Domain<escript.Domain>}         :type domain: `Domain`
1402         @param shape: suggested shape of the coefficient         :param shape: suggested shape of the coefficient
1403         @type shape: C{tuple} of C{int} values         :type shape: ``tuple`` of ``int`` values
1404         @return: the number of equations and number of solutions of the PDE if         :return: the number of equations and number of solutions of the PDE if
1405                  the coefficient has given shape. If no appropriate numbers                  the coefficient has given shape. If no appropriate numbers
1406                  could be identified, C{None} is returned                  could be identified, ``None`` is returned
1407         @rtype: C{tuple} of two C{int} values or C{None}         :rtype: ``tuple`` of two ``int`` values or ``None``
1408         """         """
1409         dim=domain.getDim()         dim=domain.getDim()
1410         if len(shape)>0:         if len(shape)>0:
# Line 320  class PDECoef(object): Line 1445  class PDECoef(object):
1445         Checks if the coefficient allows to estimate the number of solution         Checks if the coefficient allows to estimate the number of solution
1446         components.         components.
1447    
1448         @return: True if the coefficient allows an estimate of the number of         :return: True if the coefficient allows an estimate of the number of
1449                  solution components, False otherwise                  solution components, False otherwise
1450         @rtype: C{bool}         :rtype: ``bool``
1451         """         """
1452         for i in self.pattern:         for i in self.pattern:
1453               if i==self.BY_SOLUTION: return True               if i==self.BY_SOLUTION: return True
# Line 332  class PDECoef(object): Line 1457  class PDECoef(object):
1457         """         """
1458         Checks if the coefficient allows to estimate the number of equations.         Checks if the coefficient allows to estimate the number of equations.
1459    
1460         @return: True if the coefficient allows an estimate of the number of         :return: True if the coefficient allows an estimate of the number of
1461                  equations, False otherwise                  equations, False otherwise
1462         @rtype: C{bool}         :rtype: ``bool``
1463         """         """
1464         for i in self.pattern:         for i in self.pattern:
1465               if i==self.BY_EQUATION: return True               if i==self.BY_EQUATION: return True
# Line 345  class PDECoef(object): Line 1470  class PDECoef(object):
1470        Compares two tuples of possible number of equations and number of        Compares two tuples of possible number of equations and number of
1471        solutions.        solutions.
1472    
1473        @param t1: the first tuple        :param t1: the first tuple
1474        @param t2: the second tuple        :param t2: the second tuple
1475        @return: 0, 1, or -1        :return: 0, 1, or -1
1476        """        """
1477    
1478        dif=t1[0]+t1[1]-(t2[0]+t2[1])        dif=t1[0]+t1[1]-(t2[0]+t2[1])
# Line 359  class PDECoef(object): Line 1484  class PDECoef(object):
1484         """         """
1485         Builds the required shape of the coefficient.         Builds the required shape of the coefficient.
1486    
1487         @param domain: domain on which the PDE uses the coefficient         :param domain: domain on which the PDE uses the coefficient
1488         @type domain: L{Domain<escript.Domain>}         :type domain: `Domain`
1489         @param numEquations: number of equations of the PDE         :param numEquations: number of equations of the PDE
1490         @type numEquations: C{int}         :type numEquations: ``int``
1491         @param numSolutions: number of components of the PDE solution         :param numSolutions: number of components of the PDE solution
1492         @type numSolutions: C{int}         :type numSolutions: ``int``
1493         @return: shape of the coefficient         :return: shape of the coefficient
1494         @rtype: C{tuple} of C{int} values         :rtype: ``tuple`` of ``int`` values
1495         """         """
1496         dim=domain.getDim()         dim=domain.getDim()
1497         s=()         s=()
# Line 379  class PDECoef(object): Line 1504  class PDECoef(object):
1504                  s=s+(dim,)                  s=s+(dim,)
1505         return s         return s
1506    
1507    #====================================================================================================================
1508    
1509  class LinearProblem(object):  class LinearProblem(object):
1510     """     """
1511     This is the base class to define a general linear PDE-type problem for     This is the base class to define a general linear PDE-type problem for
1512     for an unknown function M{u} on a given domain defined through a     for an unknown function *u* on a given domain defined through a
1513     L{Domain<escript.Domain>} object. The problem can be given as a single     `Domain` object. The problem can be given as a single
1514     equation or as a system of equations.     equation or as a system of equations.
1515    
1516     The class assumes that some sort of assembling process is required to form     The class assumes that some sort of assembling process is required to form
1517     a problem of the form     a problem of the form
1518    
1519     M{L u=f}     *L u=f*
   
    where M{L} is an operator and M{f} is the right hand side. This operator  
    problem will be solved to get the unknown M{u}.  
   
    @cvar DEFAULT: The default method used to solve the system of linear  
                   equations  
    @cvar DIRECT: The direct solver based on LDU factorization  
    @cvar CHOLEVSKY: The direct solver based on LDLt factorization (can only be  
                     applied for symmetric PDEs)  
    @cvar PCG: The preconditioned conjugate gradient method (can only be applied  
               for symmetric PDEs)  
    @cvar CR: The conjugate residual method  
    @cvar CGS: The conjugate gradient square method  
    @cvar BICGSTAB: The stabilized BiConjugate Gradient method  
    @cvar TFQMR: Transport Free Quasi Minimal Residual method  
    @cvar MINRES: Minimum residual method  
    @cvar SSOR: The symmetric overrelaxation method  
    @cvar ILU0: The incomplete LU factorization preconditioner with no fill-in  
    @cvar ILUT: The incomplete LU factorization preconditioner with fill-in  
    @cvar JACOBI: The Jacobi preconditioner  
    @cvar GMRES: The Gram-Schmidt minimum residual method  
    @cvar PRES20: Special GMRES with restart after 20 steps and truncation after  
                  5 residuals  
    @cvar LUMPING: Matrix lumping  
    @cvar NO_REORDERING: No matrix reordering allowed  
    @cvar MINIMUM_FILL_IN: Reorder matrix to reduce fill-in during factorization  
    @cvar NESTED_DISSECTION: Reorder matrix to improve load balancing during  
                             factorization  
    @cvar PASO: PASO solver package  
    @cvar SCSL: SGI SCSL solver library  
    @cvar MKL: Intel's MKL solver library  
    @cvar UMFPACK: The UMFPACK library  
    @cvar TRILINOS: The TRILINOS parallel solver class library from Sandia Natl  
                    Labs  
    @cvar ITERATIVE: The default iterative solver  
    @cvar AMG: Algebraic Multi Grid  
    @cvar RILU: Recursive ILU  
    @cvar GS: Gauss-Seidel solver  
   
    """  
    DEFAULT= 0  
    DIRECT= 1  
    CHOLEVSKY= 2  
    PCG= 3  
    CR= 4  
    CGS= 5  
    BICGSTAB= 6  
    SSOR= 7  
    ILU0= 8  
    ILUT= 9  
    JACOBI= 10  
    GMRES= 11  
    PRES20= 12  
    LUMPING= 13  
    NO_REORDERING= 17  
    MINIMUM_FILL_IN= 18  
    NESTED_DISSECTION= 19  
    SCSL= 14  
    MKL= 15  
    UMFPACK= 16  
    ITERATIVE= 20  
    PASO= 21  
    AMG= 22  
    RILU = 23  
    TRILINOS = 24  
    NONLINEAR_GMRES = 25  
    TFQMR = 26  
    MINRES = 27  
    GS=28  
   
    SMALL_TOLERANCE=1.e-13  
    PACKAGE_KEY="package"  
    METHOD_KEY="method"  
    SYMMETRY_KEY="symmetric"  
    TOLERANCE_KEY="tolerance"  
    PRECONDITIONER_KEY="preconditioner"  
1520    
1521       where *L* is an operator and *f* is the right hand side. This operator
1522       problem will be solved to get the unknown *u*.
1523    
1524       """
1525     def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):     def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):
1526       """       """
1527       Initializes a linear problem.       Initializes a linear problem.
1528    
1529       @param domain: domain of the PDE       :param domain: domain of the PDE
1530       @type domain: L{Domain<escript.Domain>}       :type domain: `Domain`
1531       @param numEquations: number of equations. If C{None} the number of       :param numEquations: number of equations. If ``None`` the number of
1532                            equations is extracted from the coefficients.                            equations is extracted from the coefficients.
1533       @param numSolutions: number of solution components. If C{None} the number       :param numSolutions: number of solution components. If ``None`` the number
1534                            of solution components is extracted from the                            of solution components is extracted from the
1535                            coefficients.                            coefficients.
1536       @param debug: if True debug information is printed       :param debug: if True debug information is printed
1537    
1538       """       """
1539       super(LinearProblem, self).__init__()       super(LinearProblem, self).__init__()
# Line 491  class LinearProblem(object): Line 1545  class LinearProblem(object):
1545       self.__altered_coefficients=False       self.__altered_coefficients=False
1546       self.__reduce_equation_order=False       self.__reduce_equation_order=False
1547       self.__reduce_solution_order=False       self.__reduce_solution_order=False
1548       self.__tolerance=1.e-8       self.__sym=False
1549       self.__solver_method=self.DEFAULT       self.setSolverOptions()
      self.__solver_package=self.DEFAULT  
      self.__preconditioner=self.DEFAULT  
1550       self.__is_RHS_valid=False       self.__is_RHS_valid=False
1551       self.__is_operator_valid=False       self.__is_operator_valid=False
      self.__sym=False  
1552       self.__COEFFICIENTS={}       self.__COEFFICIENTS={}
1553         self.__solution_rtol=1.e99
1554         self.__solution_atol=1.e99
1555         self.setSymmetryOff()
1556       # initialize things:       # initialize things:
1557       self.resetAllCoefficients()       self.resetAllCoefficients()
1558       self.__system_type=None       self.initializeSystem()
      self.updateSystemType()  
1559     # ==========================================================================     # ==========================================================================
1560     #    general stuff:     #    general stuff:
1561     # ==========================================================================     # ==========================================================================
# Line 510  class LinearProblem(object): Line 1563  class LinearProblem(object):
1563       """       """
1564       Returns a string representation of the PDE.       Returns a string representation of the PDE.
1565    
1566       @return: a simple representation of the PDE       :return: a simple representation of the PDE
1567       @rtype: C{str}       :rtype: ``str``
1568       """       """
1569       return "<LinearProblem %d>"%id(self)       return "<LinearProblem %d>"%id(self)
1570     # ==========================================================================     # ==========================================================================
# Line 531  class LinearProblem(object): Line 1584  class LinearProblem(object):
1584    
1585     def setDebug(self, flag):     def setDebug(self, flag):
1586       """       """
1587       Switches debug output on if C{flag} is True otherwise it is switched off.       Switches debug output on if ``flag`` is True otherwise it is switched off.
1588    
1589       @param flag: desired debug status       :param flag: desired debug status
1590       @type flag: C{bool}       :type flag: ``bool``
1591       """       """
1592       if flag:       if flag:
1593           self.setDebugOn()           self.setDebugOn()
# Line 545  class LinearProblem(object): Line 1598  class LinearProblem(object):
1598       """       """
1599       Prints the text message if debug mode is switched on.       Prints the text message if debug mode is switched on.
1600    
1601       @param text: message to be printed       :param text: message to be printed
1602       @type text: C{string}       :type text: ``string``
1603       """       """
1604       if self.__debug: print "%s: %s"%(str(self),text)       if self.__debug: print "%s: %s"%(str(self),text)
1605    
# Line 557  class LinearProblem(object): Line 1610  class LinearProblem(object):
1610         """         """
1611         Introduces new coefficients into the problem.         Introduces new coefficients into the problem.
1612    
1613         Use::         Use:
1614    
1615         p.introduceCoefficients(A=PDECoef(...), B=PDECoef(...))         p.introduceCoefficients(A=PDECoef(...), B=PDECoef(...))
1616    
1617         to introduce the coefficients M{A} ans M{B}.         to introduce the coefficients *A* and *B*.
1618         """         """
1619         for name, type in coeff.items():         for name, type in coeff.items():
1620             if not isinstance(type,PDECoef):             if not isinstance(type,PDECoef):
# Line 574  class LinearProblem(object): Line 1627  class LinearProblem(object):
1627       """       """
1628       Returns the domain of the PDE.       Returns the domain of the PDE.
1629    
1630       @return: the domain of the PDE       :return: the domain of the PDE
1631       @rtype: L{Domain<escript.Domain>}       :rtype: `Domain`
1632       """       """
1633       return self.__domain       return self.__domain
1634       def getDomainStatus(self):
1635         """
1636         Return the status indicator of the domain
1637         """
1638         return self.getDomain().getStatus()
1639    
1640       def getSystemStatus(self):
1641         """
1642         Return the domain status used to build the current system
1643         """
1644         return self.__system_status
1645       def setSystemStatus(self,status=None):
1646         """
1647         Sets the system status to ``status`` if ``status`` is not present the
1648         current status of the domain is used.
1649         """
1650         if status == None:
1651             self.__system_status=self.getDomainStatus()
1652         else:
1653             self.__system_status=status
1654    
1655     def getDim(self):     def getDim(self):
1656       """       """
1657       Returns the spatial dimension of the PDE.       Returns the spatial dimension of the PDE.
1658    
1659       @return: the spatial dimension of the PDE domain       :return: the spatial dimension of the PDE domain
1660       @rtype: C{int}       :rtype: ``int``
1661       """       """
1662       return self.getDomain().getDim()       return self.getDomain().getDim()
1663    
# Line 592  class LinearProblem(object): Line 1665  class LinearProblem(object):
1665       """       """
1666       Returns the number of equations.       Returns the number of equations.
1667    
1668       @return: the number of equations       :return: the number of equations
1669       @rtype: C{int}       :rtype: ``int``
1670       @raise UndefinedPDEError: if the number of equations is not specified yet       :raise UndefinedPDEError: if the number of equations is not specified yet
1671       """       """
1672       if self.__numEquations==None:       if self.__numEquations==None:
1673           if self.__numSolutions==None:           if self.__numSolutions==None:
# Line 607  class LinearProblem(object): Line 1680  class LinearProblem(object):
1680       """       """
1681       Returns the number of unknowns.       Returns the number of unknowns.
1682    
1683       @return: the number of unknowns       :return: the number of unknowns
1684       @rtype: C{int}       :rtype: ``int``
1685       @raise UndefinedPDEError: if the number of unknowns is not specified yet       :raise UndefinedPDEError: if the number of unknowns is not specified yet
1686       """       """
1687       if self.__numSolutions==None:       if self.__numSolutions==None:
1688          if self.__numEquations==None:          if self.__numEquations==None:
# Line 622  class LinearProblem(object): Line 1695  class LinearProblem(object):
1695       """       """
1696       Returns the status of order reduction for the equation.       Returns the status of order reduction for the equation.
1697    
1698       @return: True if reduced interpolation order is used for the       :return: True if reduced interpolation order is used for the
1699                representation of the equation, False otherwise                representation of the equation, False otherwise
1700       @rtype: L{bool}       :rtype: `bool`
1701       """       """
1702       return self.__reduce_equation_order       return self.__reduce_equation_order
1703    
# Line 632  class LinearProblem(object): Line 1705  class LinearProblem(object):
1705       """       """
1706       Returns the status of order reduction for the solution.       Returns the status of order reduction for the solution.
1707    
1708       @return: True if reduced interpolation order is used for the       :return: True if reduced interpolation order is used for the
1709                representation of the solution, False otherwise                representation of the solution, False otherwise
1710       @rtype: L{bool}       :rtype: `bool`
1711       """       """
1712       return self.__reduce_solution_order       return self.__reduce_solution_order
1713    
1714     def getFunctionSpaceForEquation(self):     def getFunctionSpaceForEquation(self):
1715       """       """
1716       Returns the L{FunctionSpace<escript.FunctionSpace>} used to discretize       Returns the `FunctionSpace` used to discretize
1717       the equation.       the equation.
1718    
1719       @return: representation space of equation       :return: representation space of equation
1720       @rtype: L{FunctionSpace<escript.FunctionSpace>}       :rtype: `FunctionSpace`
1721       """       """
1722       if self.reduceEquationOrder():       if self.reduceEquationOrder():
1723           return escript.ReducedSolution(self.getDomain())           return escript.ReducedSolution(self.getDomain())
# Line 653  class LinearProblem(object): Line 1726  class LinearProblem(object):
1726    
1727     def getFunctionSpaceForSolution(self):     def getFunctionSpaceForSolution(self):
1728       """       """
1729       Returns the L{FunctionSpace<escript.FunctionSpace>} used to represent       Returns the `FunctionSpace` used to represent
1730       the solution.       the solution.
1731    
1732       @return: representation space of solution       :return: representation space of solution
1733       @rtype: L{FunctionSpace<escript.FunctionSpace>}       :rtype: `FunctionSpace`
1734       """       """
1735       if self.reduceSolutionOrder():       if self.reduceSolutionOrder():
1736           return escript.ReducedSolution(self.getDomain())           return escript.ReducedSolution(self.getDomain())
# Line 667  class LinearProblem(object): Line 1740  class LinearProblem(object):
1740     # ==========================================================================     # ==========================================================================
1741     #   solver settings:     #   solver settings:
1742     # ==========================================================================     # ==========================================================================
1743     def setSolverMethod(self,solver=None,preconditioner=None):     def setSolverOptions(self,options=None):
1744         """         """
1745         Sets a new solver method and/or preconditioner.         Sets the solver options.
   
        @param solver: new solver method to use  
        @type solver: one of L{DEFAULT}, L{ITERATIVE} L{DIRECT}, L{CHOLEVSKY},  
                      L{PCG}, L{CR}, L{CGS}, L{BICGSTAB}, L{SSOR}, L{GMRES},  
                      L{TFQMR}, L{MINRES}, L{PRES20}, L{LUMPING}, L{AMG}  
        @param preconditioner: new preconditioner to use  
        @type preconditioner: one of L{DEFAULT}, L{JACOBI} L{ILU0}, L{ILUT},  
                              L{SSOR}, L{RILU}, L{GS}  
   
        @note: the solver method chosen may not be available by the selected  
               package.  
        """  
        if solver==None: solver=self.__solver_method  
        if preconditioner==None: preconditioner=self.__preconditioner  
        if solver==None: solver=self.DEFAULT  
        if preconditioner==None: preconditioner=self.DEFAULT  
        if not (solver,preconditioner)==self.getSolverMethod():  
            self.__solver_method=solver  
            self.__preconditioner=preconditioner  
            self.updateSystemType()  
            self.trace("New solver is %s"%self.getSolverMethodName())  
   
    def getSolverMethodName(self):  
        """  
        Returns the name of the solver currently used.  
   
        @return: the name of the solver currently used  
        @rtype: C{string}  
        """  
   
        m=self.getSolverMethod()  
        p=self.getSolverPackage()  
        method=""  
        if m[0]==self.DEFAULT: method="DEFAULT"  
        elif m[0]==self.DIRECT: method= "DIRECT"  
        elif m[0]==self.ITERATIVE: method= "ITERATIVE"  
        elif m[0]==self.CHOLEVSKY: method= "CHOLEVSKY"  
        elif m[0]==self.PCG: method= "PCG"  
        elif m[0]==self.TFQMR: method= "TFQMR"  
        elif m[0]==self.MINRES: method= "MINRES"  
        elif m[0]==self.CR: method= "CR"  
        elif m[0]==self.CGS: method= "CGS"  
        elif m[0]==self.BICGSTAB: method= "BICGSTAB"  
        elif m[0]==self.SSOR: method= "SSOR"  
        elif m[0]==self.GMRES: method= "GMRES"  
        elif m[0]==self.PRES20: method= "PRES20"  
        elif m[0]==self.LUMPING: method= "LUMPING"  
        elif m[0]==self.AMG: method= "AMG"  
        if m[1]==self.DEFAULT: method+="+DEFAULT"  
        elif m[1]==self.JACOBI: method+= "+JACOBI"  
        elif m[1]==self.ILU0: method+= "+ILU0"  
        elif m[1]==self.ILUT: method+= "+ILUT"  
        elif m[1]==self.SSOR: method+= "+SSOR"  
        elif m[1]==self.AMG: method+= "+AMG"  
        elif m[1]==self.RILU: method+= "+RILU"  
        elif m[1]==self.GS: method+= "+GS"  
        if p==self.DEFAULT: package="DEFAULT"  
        elif p==self.PASO: package= "PASO"  
        elif p==self.MKL: package= "MKL"  
        elif p==self.SCSL: package= "SCSL"  
        elif p==self.UMFPACK: package= "UMFPACK"  
        elif p==self.TRILINOS: package= "TRILINOS"  
        else : method="unknown"  
        return "%s solver of %s package"%(method,package)  
   
    def getSolverMethod(self):  
        """  
        Returns the solver method and preconditioner used.  
   
        @return: the solver method currently used.  
        @rtype: C{tuple} of C{int}  
        """  
        return self.__solver_method,self.__preconditioner  
   
    def setSolverPackage(self,package=None):  
        """  
        Sets a new solver package.  
   
        @param package: the new solver package  
        @type package: one of L{DEFAULT}, L{PASO} L{SCSL}, L{MKL}, L{UMFPACK},  
                       L{TRILINOS}  
        """  
        if package==None: package=self.DEFAULT  
        if not package==self.getSolverPackage():  
            self.__solver_package=package  
            self.updateSystemType()  
            self.trace("New solver is %s"%self.getSolverMethodName())  
   
    def getSolverPackage(self):  
        """  
        Returns the package of the solver.  
   
        @return: the solver package currently being used  
        @rtype: C{int}  
        """  
        return self.__solver_package  
1746    
1747           :param options: the new solver options. If equal ``None``, the solver options are set to the default.
1748           :type options: `SolverOptions` or ``None``
1749           :note: The symmetry flag of options is overwritten by the symmetry flag of the `LinearProblem`.
1750           """
1751           if options==None:
1752              self.__solver_options=SolverOptions()
1753           elif isinstance(options, SolverOptions):
1754              self.__solver_options=options
1755           else:
1756              raise ValueError,"options must be a SolverOptions object."
1757           self.__solver_options.setSymmetry(self.__sym)
1758        
1759       def getSolverOptions(self):
1760           """
1761           Returns the solver options
1762      
1763           :rtype: `SolverOptions`
1764           """
1765           self.__solver_options.setSymmetry(self.__sym)
1766           return self.__solver_options
1767          
1768     def isUsingLumping(self):     def isUsingLumping(self):
1769        """        """
1770        Checks if matrix lumping is the current solver method.        Checks if matrix lumping is the current solver method.
1771    
1772        @return: True if the current solver method is lumping        :return: True if the current solver method is lumping
1773        @rtype: C{bool}        :rtype: ``bool``
1774        """        """
1775        return self.getSolverMethod()[0]==self.LUMPING        return self.getSolverOptions().getSolverMethod() in [ SolverOptions.ROWSUM_LUMPING, SolverOptions.HRZ_LUMPING ]
   
    def setTolerance(self,tol=1.e-8):  
        """  
        Resets the tolerance for the solver method to C{tol}.  
   
        @param tol: new tolerance for the solver. If C{tol} is lower than the  
                    current tolerence the system will be resolved.  
        @type tol: positive C{float}  
        @raise ValueError: if tolerance is not positive  
        """  
        if not tol>0:  
            raise ValueError,"Tolerance has to be positive"  
        if tol<self.getTolerance(): self.invalidateSolution()  
        self.trace("New tolerance %e"%tol)  
        self.__tolerance=tol  
        return  
   
    def getTolerance(self):  
        """  
        Returns the tolerance currently set for the solution.  
   
        @return: tolerance currently used  
        @rtype: C{float}  
        """  
        return self.__tolerance  
   
1776     # ==========================================================================     # ==========================================================================
1777     #    symmetry  flag:     #    symmetry  flag:
1778     # ==========================================================================     # ==========================================================================
# Line 808  class LinearProblem(object): Line 1780  class LinearProblem(object):
1780        """        """
1781        Checks if symmetry is indicated.        Checks if symmetry is indicated.
1782    
1783        @return: True if a symmetric PDE is indicated, False otherwise        :return: True if a symmetric PDE is indicated, False otherwise
1784        @rtype: C{bool}        :rtype: ``bool``
1785          :note: the method is equivalent to use getSolverOptions().isSymmetric()
1786        """        """
1787        return self.__sym        self.getSolverOptions().isSymmetric()
1788    
1789     def setSymmetryOn(self):     def setSymmetryOn(self):
1790        """        """
1791        Sets the symmetry flag.        Sets the symmetry flag.
1792          :note: The method overwrites the symmetry flag set by the solver options
1793        """        """
1794        if not self.isSymmetric():        self.__sym=True
1795           self.trace("PDE is set to be symmetric")        self.getSolverOptions().setSymmetryOn()
          self.__sym=True  
          self.updateSystemType()  
1796    
1797     def setSymmetryOff(self):     def setSymmetryOff(self):
1798        """        """
1799        Clears the symmetry flag.        Clears the symmetry flag.
1800          :note: The method overwrites the symmetry flag set by the solver options
1801        """        """
1802        if self.isSymmetric():        self.__sym=False
1803           self.trace("PDE is set to be nonsymmetric")        self.getSolverOptions().setSymmetryOff()
          self.__sym=False  
          self.updateSystemType()  
1804    
1805     def setSymmetryTo(self,flag=False):     def setSymmetry(self,flag=False):
1806        """        """
1807        Sets the symmetry flag to C{flag}.        Sets the symmetry flag to ``flag``.
1808    
1809        @param flag: If True, the symmetry flag is set otherwise reset.        :param flag: If True, the symmetry flag is set otherwise reset.
1810        @type flag: C{bool}        :type flag: ``bool``
1811          :note: The method overwrites the symmetry flag set by the solver options
1812        """        """
1813        if flag:        self.getSolverOptions().setSymmetry(flag)
          self.setSymmetryOn()  
       else:  
          self.setSymmetryOff()  
   
1814     # ==========================================================================     # ==========================================================================
1815     # function space handling for the equation as well as the solution     # function space handling for the equation as well as the solution
1816     # ==========================================================================     # ==========================================================================
# Line 850  class LinearProblem(object): Line 1818  class LinearProblem(object):
1818       """       """
1819       Switches reduced order on for solution and equation representation.       Switches reduced order on for solution and equation representation.
1820    
1821       @raise RuntimeError: if order reduction is altered after a coefficient has       :raise RuntimeError: if order reduction is altered after a coefficient has
1822                            been set                            been set
1823       """       """
1824       self.setReducedOrderForSolutionOn()       self.setReducedOrderForSolutionOn()
# Line 860  class LinearProblem(object): Line 1828  class LinearProblem(object):
1828       """       """
1829       Switches reduced order off for solution and equation representation       Switches reduced order off for solution and equation representation
1830    
1831       @raise RuntimeError: if order reduction is altered after a coefficient has       :raise RuntimeError: if order reduction is altered after a coefficient has
1832                            been set                            been set
1833       """       """
1834       self.setReducedOrderForSolutionOff()       self.setReducedOrderForSolutionOff()
# Line 871  class LinearProblem(object): Line 1839  class LinearProblem(object):
1839       Sets order reduction state for both solution and equation representation       Sets order reduction state for both solution and equation representation
1840       according to flag.       according to flag.
1841    
1842       @param flag: if True, the order reduction is switched on for both solution       :param flag: if True, the order reduction is switched on for both solution
1843                    and equation representation, otherwise or if flag is not                    and equation representation, otherwise or if flag is not
1844                    present order reduction is switched off                    present order reduction is switched off
1845       @type flag: C{bool}       :type flag: ``bool``
1846       @raise RuntimeError: if order reduction is altered after a coefficient has       :raise RuntimeError: if order reduction is altered after a coefficient has
1847                            been set                            been set
1848       """       """
1849       self.setReducedOrderForSolutionTo(flag)       self.setReducedOrderForSolutionTo(flag)
# Line 886  class LinearProblem(object): Line 1854  class LinearProblem(object):
1854       """       """
1855       Switches reduced order on for solution representation.       Switches reduced order on for solution representation.
1856    
1857       @raise RuntimeError: if order reduction is altered after a coefficient has       :raise RuntimeError: if order reduction is altered after a coefficient has
1858                            been set                            been set
1859       """       """
1860       if not self.__reduce_solution_order:       if not self.__reduce_solution_order:
# Line 900  class LinearProblem(object): Line 1868  class LinearProblem(object):
1868       """       """
1869       Switches reduced order off for solution representation       Switches reduced order off for solution representation
1870    
1871       @raise RuntimeError: if order reduction is altered after a coefficient has       :raise RuntimeError: if order reduction is altered after a coefficient has
1872                            been set.                            been set.
1873       """       """
1874       if self.__reduce_solution_order:       if self.__reduce_solution_order:
# Line 914  class LinearProblem(object): Line 1882  class LinearProblem(object):
1882       """       """
1883       Sets order reduction state for solution representation according to flag.       Sets order reduction state for solution representation according to flag.
1884    
1885       @param flag: if flag is True, the order reduction is switched on for       :param flag: if flag is True, the order reduction is switched on for
1886                    solution representation, otherwise or if flag is not present                    solution representation, otherwise or if flag is not present
1887                    order reduction is switched off                    order reduction is switched off
1888       @type flag: C{bool}       :type flag: ``bool``
1889       @raise RuntimeError: if order reduction is altered after a coefficient has       :raise RuntimeError: if order reduction is altered after a coefficient has
1890                            been set                            been set
1891       """       """
1892       if flag:       if flag:
# Line 930  class LinearProblem(object): Line 1898  class LinearProblem(object):
1898       """       """
1899       Switches reduced order on for equation representation.       Switches reduced order on for equation representation.
1900    
1901       @raise RuntimeError: if order reduction is altered after a coefficient has       :raise RuntimeError: if order reduction is altered after a coefficient has
1902                            been set                            been set
1903       """       """
1904       if not self.__reduce_equation_order:       if not self.__reduce_equation_order:
# Line 944  class LinearProblem(object): Line 1912  class LinearProblem(object):
1912       """       """
1913       Switches reduced order off for equation representation.       Switches reduced order off for equation representation.
1914    
1915       @raise RuntimeError: if order reduction is altered after a coefficient has       :raise RuntimeError: if order reduction is altered after a coefficient has
1916                            been set                            been set
1917       """       """
1918       if self.__reduce_equation_order:       if self.__reduce_equation_order:
# Line 958  class LinearProblem(object): Line 1926  class LinearProblem(object):
1926       """       """
1927       Sets order reduction state for equation representation according to flag.       Sets order reduction state for equation representation according to flag.
1928    
1929       @param flag: if flag is True, the order reduction is switched on for       :param flag: if flag is True, the order reduction is switched on for
1930                    equation representation, otherwise or if flag is not present                    equation representation, otherwise or if flag is not present
1931                    order reduction is switched off                    order reduction is switched off
1932       @type flag: C{bool}       :type flag: ``bool``
1933       @raise RuntimeError: if order reduction is altered after a coefficient has       :raise RuntimeError: if order reduction is altered after a coefficient has
1934                            been set                            been set
1935       """       """
1936       if flag:       if flag:
1937          self.setReducedOrderForEquationOn()          self.setReducedOrderForEquationOn()
1938       else:       else:
1939          self.setReducedOrderForEquationOff()          self.setReducedOrderForEquationOff()
1940       def getOperatorType(self):
    def updateSystemType(self):  
      """  
      Reassesses the system type and, if a new matrix is needed, resets the  
      system.  
      """  
      new_system_type=self.getRequiredSystemType()  
      if not new_system_type==self.__system_type:  
          self.trace("Matrix type is now %d."%new_system_type)  
          self.__system_type=new_system_type  
          self.initializeSystem()  
   
    def getSystemType(self):  
1941        """        """
1942        Returns the current system type.        Returns the current system type.
1943        """        """
1944        return self.__system_type        return self.__operator_type
1945    
1946     def checkSymmetricTensor(self,name,verbose=True):     def checkSymmetricTensor(self,name,verbose=True):
1947        """        """
1948        Tests a coefficient for symmetry.        Tests a coefficient for symmetry.
1949    
1950        @param name: name of the coefficient        :param name: name of the coefficient
1951        @type name: C{str}        :type name: ``str``
1952        @param verbose: if set to True or not present a report on coefficients        :param verbose: if set to True or not present a report on coefficients
1953                        which break the symmetry is printed.                        which break the symmetry is printed.
1954        @type verbose: C{bool}        :type verbose: ``bool``
1955        @return: True if coefficient C{name} is symmetric        :return: True if coefficient ``name`` is symmetric
1956        @rtype: C{bool}        :rtype: ``bool``
1957        """        """
1958          SMALL_TOLERANCE=util.EPSILON*10.
1959        A=self.getCoefficient(name)        A=self.getCoefficient(name)
1960        verbose=verbose or self.__debug        verbose=verbose or self.__debug
1961        out=True        out=True
1962        if not A.isEmpty():        if not A.isEmpty():
1963           tol=util.Lsup(A)*self.SMALL_TOLERANCE           tol=util.Lsup(A)*SMALL_TOLERANCE
1964           s=A.getShape()           s=A.getShape()
1965           if A.getRank() == 4:           if A.getRank() == 4:
1966              if s[0]==s[2] and s[1] == s[3]:              if s[0]==s[2] and s[1] == s[3]:
# Line 1037  class LinearProblem(object): Line 1994  class LinearProblem(object):
1994        """        """
1995        Tests two coefficients for reciprocal symmetry.        Tests two coefficients for reciprocal symmetry.
1996    
1997        @param name0: name of the first coefficient        :param name0: name of the first coefficient
1998        @type name0: C{str}        :type name0: ``str``
1999        @param name1: name of the second coefficient        :param name1: name of the second coefficient
2000        @type name1: C{str}        :type name1: ``str``
2001        @param verbose: if set to True or not present a report on coefficients        :param verbose: if set to True or not present a report on coefficients
2002                        which break the symmetry is printed                        which break the symmetry is printed
2003        @type verbose: C{bool}        :type verbose: ``bool``
2004        @return: True if coefficients C{name0} and C{name1} are reciprocally        :return: True if coefficients ``name0`` and ``name1`` are reciprocally
2005                 symmetric.                 symmetric.
2006        @rtype: C{bool}        :rtype: ``bool``
2007        """        """
2008          SMALL_TOLERANCE=util.EPSILON*10.
2009        B=self.getCoefficient(name0)        B=self.getCoefficient(name0)
2010        C=self.getCoefficient(name1)        C=self.getCoefficient(name1)
2011        verbose=verbose or self.__debug        verbose=verbose or self.__debug
# Line 1061  class LinearProblem(object): Line 2019  class LinearProblem(object):
2019        elif not B.isEmpty() and not C.isEmpty():        elif not B.isEmpty() and not C.isEmpty():
2020           sB=B.getShape()           sB=B.getShape()
2021           sC=C.getShape()           sC=C.getShape()
2022           tol=(util.Lsup(B)+util.Lsup(C))*self.SMALL_TOLERANCE/2.           tol=(util.Lsup(B)+util.Lsup(C))*SMALL_TOLERANCE/2.
2023           if len(sB) != len(sC):           if len(sB) != len(sC):
2024               if verbose: print "non-symmetric problem because ranks of %s (=%s) and %s (=%s) are different."%(name0,len(sB),name1,len(sC))               if verbose: print "non-symmetric problem because ranks of %s (=%s) and %s (=%s) are different."%(name0,len(sB),name1,len(sC))
2025               out=False               out=False
# Line 1094  class LinearProblem(object): Line 2052  class LinearProblem(object):
2052    
2053     def getCoefficient(self,name):     def getCoefficient(self,name):
2054       """       """
2055       Returns the value of the coefficient C{name}.       Returns the value of the coefficient ``name``.
2056    
2057       @param name: name of the coefficient requested       :param name: name of the coefficient requested
2058       @type name: C{string}       :type name: ``string``
2059       @return: the value of the coefficient       :return: the value of the coefficient
2060       @rtype: L{Data<escript.Data>}       :rtype: `Data`
2061       @raise IllegalCoefficient: if C{name} is not a coefficient of the PDE       :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE
2062       """       """
2063       if self.hasCoefficient(name):       if self.hasCoefficient(name):
2064           return self.__COEFFICIENTS[name].getValue()           return self.__COEFFICIENTS[name].getValue()
# Line 1109  class LinearProblem(object): Line 2067  class LinearProblem(object):
2067    
2068     def hasCoefficient(self,name):     def hasCoefficient(self,name):
2069       """       """
2070       Returns True if C{name} is the name of a coefficient.       Returns True if ``name`` is the name of a coefficient.
2071    
2072       @param name: name of the coefficient enquired       :param name: name of the coefficient enquired
2073       @type name: C{string}       :type name: ``string``
2074       @return: True if C{name} is the name of a coefficient of the general PDE,       :return: True if ``name`` is the name of a coefficient of the general PDE,
2075                False otherwise                False otherwise
2076       @rtype: C{bool}       :rtype: ``bool``
2077       """       """
2078       return self.__COEFFICIENTS.has_key(name)       return self.__COEFFICIENTS.has_key(name)
2079    
2080     def createCoefficient(self, name):     def createCoefficient(self, name):
2081       """       """
2082       Creates a L{Data<escript.Data>} object corresponding to coefficient       Creates a `Data` object corresponding to coefficient
2083       C{name}.       ``name``.
2084    
2085       @return: the coefficient C{name} initialized to 0       :return: the coefficient ``name`` initialized to 0
2086       @rtype: L{Data<escript.Data>}       :rtype: `Data`
2087       @raise IllegalCoefficient: if C{name} is not a coefficient of the PDE       :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE
2088       """       """
2089       if self.hasCoefficient(name):       if self.hasCoefficient(name):
2090          return escript.Data(0.,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name))          return escript.Data(0.,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name))
# Line 1135  class LinearProblem(object): Line 2093  class LinearProblem(object):
2093    
2094     def getFunctionSpaceForCoefficient(self,name):     def getFunctionSpaceForCoefficient(self,name):
2095       """       """
2096       Returns the L{FunctionSpace<escript.FunctionSpace>} to be used for       Returns the `FunctionSpace` to be used for
2097       coefficient C{name}.       coefficient ``name``.
2098    
2099       @param name: name of the coefficient enquired       :param name: name of the coefficient enquired
2100       @type name: C{string}       :type name: ``string``
2101       @return: the function space to be used for coefficient C{name}       :return: the function space to be used for coefficient ``name``
2102       @rtype: L{FunctionSpace<escript.FunctionSpace>}       :rtype: `FunctionSpace`
2103       @raise IllegalCoefficient: if C{name} is not a coefficient of the PDE       :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE
2104       """       """
2105       if self.hasCoefficient(name):       if self.hasCoefficient(name):
2106          return self.__COEFFICIENTS[name].getFunctionSpace(self.getDomain())          return self.__COEFFICIENTS[name].getFunctionSpace(self.getDomain())
# Line 1151  class LinearProblem(object): Line 2109  class LinearProblem(object):
2109    
2110     def getShapeOfCoefficient(self,name):     def getShapeOfCoefficient(self,name):
2111       """       """
2112       Returns the shape of the coefficient C{name}.       Returns the shape of the coefficient ``name``.
2113    
2114       @param name: name of the coefficient enquired       :param name: name of the coefficient enquired
2115       @type name: C{string}       :type name: ``string``
2116       @return: the shape of the coefficient C{name}       :return: the shape of the coefficient ``name``
2117       @rtype: C{tuple} of C{int}       :rtype: ``tuple`` of ``int``
2118       @raise IllegalCoefficient: if C{name} is not a coefficient of the PDE       :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE
2119       """       """
2120       if self.hasCoefficient(name):       if self.hasCoefficient(name):
2121          return self.__COEFFICIENTS[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())          return self.__COEFFICIENTS[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())
# Line 1173  class LinearProblem(object): Line 2131  class LinearProblem(object):
2131    
2132     def alteredCoefficient(self,name):     def alteredCoefficient(self,name):
2133       """       """
2134       Announces that coefficient C{name} has been changed.       Announces that coefficient ``name`` has been changed.
2135    
2136       @param name: name of the coefficient affected       :param name: name of the coefficient affected
2137       @type name: C{string}       :type name: ``string``
2138       @raise IllegalCoefficient: if C{name} is not a coefficient of the PDE       :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE
2139       @note: if C{name} is q or r, the method will not trigger a rebuild of the       :note: if ``name`` is q or r, the method will not trigger a rebuild of the
2140              system as constraints are applied to the solved system.              system as constraints are applied to the solved system.
2141       """       """
2142       if self.hasCoefficient(name):       if self.hasCoefficient(name):
# Line 1206  class LinearProblem(object): Line 2164  class LinearProblem(object):
2164         """         """
2165         Returns True if the solution is still valid.         Returns True if the solution is still valid.
2166         """         """
2167           if not self.getDomainStatus()==self.getSystemStatus(): self.invalidateSolution()
2168           if self.__solution_rtol>self.getSolverOptions().getTolerance() or \
2169              self.__solution_atol>self.getSolverOptions().getAbsoluteTolerance():
2170             self.invalidateSolution()  
2171         return self.__is_solution_valid         return self.__is_solution_valid
2172    
2173     def validOperator(self):     def validOperator(self):
# Line 1226  class LinearProblem(object): Line 2188  class LinearProblem(object):
2188         """         """
2189         Returns True if the operator is still valid.         Returns True if the operator is still valid.
2190         """         """
2191           if not self.getDomainStatus()==self.getSystemStatus(): self.invalidateOperator()
2192           if not self.getRequiredOperatorType()==self.getOperatorType(): self.invalidateOperator()
2193         return self.__is_operator_valid         return self.__is_operator_valid
2194    
2195     def validRightHandSide(self):     def validRightHandSide(self):
# Line 1238  class LinearProblem(object): Line 2202  class LinearProblem(object):
2202         """         """
2203         Indicates the right hand side has to be rebuilt next time it is used.         Indicates the right hand side has to be rebuilt next time it is used.
2204         """         """
2205         if self.isRightHandSideValid(): self.trace("Right hand side has to be rebuilt.")         self.trace("Right hand side has to be rebuilt.")
2206         self.invalidateSolution()         self.invalidateSolution()
2207         self.__is_RHS_valid=False         self.__is_RHS_valid=False
2208    
# Line 1246  class LinearProblem(object): Line 2210  class LinearProblem(object):
2210         """         """
2211         Returns True if the operator is still valid.         Returns True if the operator is still valid.
2212         """         """
2213           if not self.getDomainStatus()==self.getSystemStatus(): self.invalidateRightHandSide()
2214         return self.__is_RHS_valid         return self.__is_RHS_valid
2215    
2216     def invalidateSystem(self):     def invalidateSystem(self):
2217         """         """
2218         Announces that everything has to be rebuilt.         Announces that everything has to be rebuilt.
2219         """         """
        if self.isRightHandSideValid(): self.trace("System has to be rebuilt.")  
2220         self.invalidateSolution()         self.invalidateSolution()
2221         self.invalidateOperator()         self.invalidateOperator()
2222         self.invalidateRightHandSide()         self.invalidateRightHandSide()
# Line 1268  class LinearProblem(object): Line 2232  class LinearProblem(object):
2232         Resets the system clearing the operator, right hand side and solution.         Resets the system clearing the operator, right hand side and solution.
2233         """         """
2234         self.trace("New System has been created.")         self.trace("New System has been created.")
2235           self.__operator_type=None
2236           self.setSystemStatus()
2237         self.__operator=escript.Operator()         self.__operator=escript.Operator()
2238         self.__righthandside=escript.Data()         self.__righthandside=escript.Data()
2239         self.__solution=escript.Data()         self.__solution=escript.Data()
# Line 1277  class LinearProblem(object): Line 2243  class LinearProblem(object):
2243       """       """
2244       Returns the operator of the linear problem.       Returns the operator of the linear problem.
2245    
2246       @return: the operator of the problem       :return: the operator of the problem
2247       """       """
2248       return self.getSystem()[0]       return self.getSystem()[0]
2249    
# Line 1285  class LinearProblem(object): Line 2251  class LinearProblem(object):
2251       """       """
2252       Returns the right hand side of the linear problem.       Returns the right hand side of the linear problem.
2253    
2254       @return: the right hand side of the problem       :return: the right hand side of the problem
2255       @rtype: L{Data<escript.Data>}       :rtype: `Data`
2256       """       """
2257       return self.getSystem()[1]       return self.getSystem()[1]
2258    
# Line 1320  class LinearProblem(object): Line 2286  class LinearProblem(object):
2286             self.__solution.setToZero()             self.__solution.setToZero()
2287             self.trace("Solution is reset to zero.")             self.trace("Solution is reset to zero.")
2288    
2289     def setSolution(self,u):     def setSolution(self,u, validate=True):
2290         """         """
2291         Sets the solution assuming that makes the system valid.         Sets the solution assuming that makes the system valid with the tolrance
2292           defined by the solver options
2293         """         """
2294           if validate:
2295          self.__solution_rtol=self.getSolverOptions().getTolerance()
2296          self.__solution_atol=self.getSolverOptions().getAbsoluteTolerance()
2297          self.validSolution()
2298         self.__solution=u         self.__solution=u
        self.validSolution()  
   
2299     def getCurrentSolution(self):     def getCurrentSolution(self):
2300         """         """
2301         Returns the solution in its current state.         Returns the solution in its current state.
# Line 1356  class LinearProblem(object): Line 2325  class LinearProblem(object):
2325         Makes sure that the operator is instantiated and returns it initialized         Makes sure that the operator is instantiated and returns it initialized
2326         with zeros.         with zeros.
2327         """         """
2328         if self.__operator.isEmpty():         if self.getOperatorType() == None:
2329             if self.isUsingLumping():             if self.isUsingLumping():
2330                 self.__operator=self.createSolution()                 self.__operator=self.createSolution()
2331             else:             else:
2332                 self.__operator=self.createOperator()                 self.__operator=self.createOperator()
2333           self.__operator_type=self.getRequiredOperatorType()
2334         else:         else:
2335             if self.isUsingLumping():             if self.isUsingLumping():
2336                 self.__operator.setToZero()                 self.__operator.setToZero()
2337             else:             else:
2338                 self.__operator.resetValues()                 if self.getOperatorType() == self.getRequiredOperatorType():
2339                       self.__operator.resetValues()
2340                   else:
2341                       self.__operator=self.createOperator()
2342                   self.__operator_type=self.getRequiredOperatorType()
2343             self.trace("Operator reset to zero")             self.trace("Operator reset to zero")
2344    
2345     def getCurrentOperator(self):     def getCurrentOperator(self):
2346         """         """
2347         Returns the operator in its current state.         Returns the operator in its current state.
2348         """         """
        if self.__operator.isEmpty(): self.resetOperator()  
2349         return self.__operator         return self.__operator
2350    
2351     def setValue(self,**coefficients):     def setValue(self,**coefficients):
2352        """        """
2353        Sets new values to coefficients.        Sets new values to coefficients.
2354    
2355        @raise IllegalCoefficient: if an unknown coefficient keyword is used        :raise IllegalCoefficient: if an unknown coefficient keyword is used
2356        """        """
2357        # check if the coefficients are  legal:        # check if the coefficients are  legal:
2358        for i in coefficients.iterkeys():        for i in coefficients.iterkeys():
# Line 1393  class LinearProblem(object): Line 2366  class LinearProblem(object):
2366              elif hasattr(d,"getShape"):              elif hasattr(d,"getShape"):
2367                  s=d.getShape()                  s=d.getShape()
2368              else:              else:
2369                  s=numarray.array(d).shape                  s=numpy.array(d).shape
2370              if s!=None:              if s!=None:
2371                  # get number of equations and number of unknowns:                  # get number of equations and number of unknowns:
2372                  res=self.__COEFFICIENTS[i].estimateNumEquationsAndNumSolutions(self.getDomain(),s)                  res=self.__COEFFICIENTS[i].estimateNumEquationsAndNumSolutions(self.getDomain(),s)
# Line 1434  class LinearProblem(object): Line 2407  class LinearProblem(object):
2407     # methods that are typically overwritten when implementing a particular     # methods that are typically overwritten when implementing a particular
2408     # linear problem     # linear problem
2409     # ==========================================================================     # ==========================================================================
2410     def getRequiredSystemType(self):     def getRequiredOperatorType(self):
2411        """        """
2412        Returns the system type which needs to be used by the current set up.        Returns the system type which needs to be used by the current set up.
2413    
2414        @note: Typically this method is overwritten when implementing a        :note: Typically this method is overwritten when implementing a
2415               particular linear problem.               particular linear problem.
2416        """        """
2417        return 0        return None
2418    
2419     def createOperator(self):     def createOperator(self):
2420         """         """
2421         Returns an instance of a new operator.         Returns an instance of a new operator.
2422    
2423         @note: This method is overwritten when implementing a particular         :note: This method is overwritten when implementing a particular
2424                linear problem.                linear problem.
2425         """         """
2426         return escript.Operator()         return escript.Operator()
# Line 1456  class LinearProblem(object): Line 2429  class LinearProblem(object):
2429        """        """
2430        Tests the PDE for symmetry.        Tests the PDE for symmetry.
2431    
2432        @param verbose: if set to True or not present a report on coefficients        :param verbose: if set to True or not present a report on coefficients
2433                        which break the symmetry is printed                        which break the symmetry is printed
2434        @type verbose: C{bool}        :type verbose: ``bool``
2435        @return: True if the problem is symmetric        :return: True if the problem is symmetric
2436        @rtype: C{bool}        :rtype: ``bool``
2437        @note: Typically this method is overwritten when implementing a        :note: Typically this method is overwritten when implementing a
2438               particular linear problem.               particular linear problem.
2439        """        """
2440        out=True        out=True
# Line 1471  class LinearProblem(object): Line 2444  class LinearProblem(object):
2444         """         """
2445         Returns the solution of the problem.         Returns the solution of the problem.
2446    
2447         @return: the solution         :return: the solution
2448         @rtype: L{Data<escript.Data>}         :rtype: `Data`
2449    
2450         @note: This method is overwritten when implementing a particular         :note: This method is overwritten when implementing a particular
2451                linear problem.                linear problem.
2452         """         """
2453         return self.getCurrentSolution()         return self.getCurrentSolution()
# Line 1483  class LinearProblem(object): Line 2456  class LinearProblem(object):
2456         """         """
2457         Returns the operator and right hand side of the PDE.         Returns the operator and right hand side of the PDE.
2458    
2459         @return: the discrete version of the PDE         :return: the discrete version of the PDE
2460         @rtype: C{tuple} of L{Operator,<escript.Operator>} and L{Data<escript.Data>}.         :rtype: ``tuple`` of `Operator` and `Data`.
2461    
2462         @note: This method is overwritten when implementing a particular         :note: This method is overwritten when implementing a particular
2463                linear problem.                linear problem.
2464         """         """
2465         return (self.getCurrentOperator(), self.getCurrentRightHandSide())         return (self.getCurrentOperator(), self.getCurrentRightHandSide())
 #=====================  
2466    
2467  class LinearPDE(LinearProblem):  class LinearPDE(LinearProblem):
2468     """     """
2469     This class is used to define a general linear, steady, second order PDE     This class is used to define a general linear, steady, second order PDE
2470     for an unknown function M{u} on a given domain defined through a     for an unknown function *u* on a given domain defined through a
2471     L{Domain<escript.Domain>} object.     `Domain` object.
2472    
2473     For a single PDE having a solution with a single component the linear PDE     For a single PDE having a solution with a single component the linear PDE
2474     is defined in the following form:     is defined in the following form:
2475    
2476     M{-(grad(A[j,l]+A_reduced[j,l])*grad(u)[l]+(B[j]+B_reduced[j])u)[j]+(C[l]+C_reduced[l])*grad(u)[l]+(D+D_reduced)=-grad(X+X_reduced)[j,j]+(Y+Y_reduced)}     *-(grad(A[j,l]+A_reduced[j,l])*grad(u)[l]+(B[j]+B_reduced[j])u)[j]+(C[l]+C_reduced[l])*grad(u)[l]+(D+D_reduced)=-grad(X+X_reduced)[j,j]+(Y+Y_reduced)*
2477    
2478     where M{grad(F)} denotes the spatial derivative of M{F}. Einstein's     where *grad(F)* denotes the spatial derivative of *F*. Einstein's
2479     summation convention, ie. summation over indexes appearing twice in a term     summation convention, ie. summation over indexes appearing twice in a term
2480     of a sum performed, is used.     of a sum performed, is used.
2481     The coefficients M{A}, M{B}, M{C}, M{D}, M{X} and M{Y} have to be specified     The coefficients *A*, *B*, *C*, *D*, *X* and *Y* have to be specified
2482     through L{Data<escript.Data>} objects in L{Function<escript.Function>} and     through `Data` objects in `Function` and
2483     the coefficients M{A_reduced}, M{B_reduced}, M{C_reduced}, M{D_reduced},     the coefficients *A_reduced*, *B_reduced*, *C_reduced*, *D_reduced*,
2484     M{X_reduced} and M{Y_reduced} have to be specified through     *X_reduced* and *Y_reduced* have to be specified through
2485     L{Data<escript.Data>} objects in L{ReducedFunction<escript.ReducedFunction>}.     `Data` objects in `ReducedFunction`.
2486     It is also allowed to use objects that can be converted into such     It is also allowed to use objects that can be converted into such
2487     L{Data<escript.Data>} objects. M{A} and M{A_reduced} are rank two, M{B},     `Data` objects. *A* and *A_reduced* are rank two, *B*,
2488     M{C}, M{X}, M{B_reduced}, M{C_reduced} and M{X_reduced} are rank one and     *C*, *X*, *B_reduced*, *C_reduced* and *X_reduced* are rank one and
2489     M{D}, M{D_reduced}, M{Y} and M{Y_reduced} are scalar.     *D*, *D_reduced*, *Y* and *Y_reduced* are scalar.
2490    
2491     The following natural boundary conditions are considered:     The following natural boundary conditions are considered:
2492    
2493     M{n[j]*((A[i,j]+A_reduced[i,j])*grad(u)[l]+(B+B_reduced)[j]*u)+(d+d_reduced)*u=n[j]*(X[j]+X_reduced[j])+y}     *n[j]*((A[i,j]+A_reduced[i,j])*grad(u)[l]+(B+B_reduced)[j]*u)+(d+d_reduced)*u=n[j]*(X[j]+X_reduced[j])+y*
2494    
2495     where M{n} is the outer normal field. Notice that the coefficients M{A},     where *n* is the outer normal field. Notice that the coefficients *A*,
2496     M{A_reduced}, M{B}, M{B_reduced}, M{X} and M{X_reduced} are defined in the     *A_reduced*, *B*, *B_reduced*, *X* and *X_reduced* are defined in the
2497     PDE. The coefficients M{d} and M{y} are each a scalar in     PDE. The coefficients *d* and *y* are each a scalar in
2498     L{FunctionOnBoundary<escript.FunctionOnBoundary>} and the coefficients     `FunctionOnBoundary` and the coefficients
2499     M{d_reduced} and M{y_reduced} are each a scalar in     *d_reduced* and *y_reduced* are each a scalar in
2500     L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.     `ReducedFunctionOnBoundary`.
2501    
2502     Constraints for the solution prescribe the value of the solution at certain     Constraints for the solution prescribe the value of the solution at certain
2503     locations in the domain. They have the form     locations in the domain. They have the form
2504    
2505     M{u=r} where M{q>0}     *u=r* where *q>0*
2506    
2507     M{r} and M{q} are each scalar where M{q} is the characteristic function     *r* and *q* are each scalar where *q* is the characteristic function
2508     defining where the constraint is applied. The constraints override any     defining where the constraint is applied. The constraints override any
2509     other condition set by the PDE or the boundary condition.     other condition set by the PDE or the boundary condition.
2510    
2511     The PDE is symmetrical if     The PDE is symmetrical if
2512    
2513     M{A[i,j]=A[j,i]}  and M{B[j]=C[j]} and M{A_reduced[i,j]=A_reduced[j,i]}     *A[i,j]=A[j,i]*  and *B[j]=C[j]* and *A_reduced[i,j]=A_reduced[j,i]*
2514     and M{B_reduced[j]=C_reduced[j]}     and *B_reduced[j]=C_reduced[j]*
2515    
2516     For a system of PDEs and a solution with several components the PDE has the     For a system of PDEs and a solution with several components the PDE has the
2517     form     form
2518    
2519     M{-grad((A[i,j,k,l]+A_reduced[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k])[j]+(C[i,k,l]+C_reduced[i,k,l])*grad(u[k])[l]+(D[i,k]+D_reduced[i,k]*u[k] =-grad(X[i,j]+X_reduced[i,j])[j]+Y[i]+Y_reduced[i] }     *-grad((A[i,j,k,l]+A_reduced[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k])[j]+(C[i,k,l]+C_reduced[i,k,l])*grad(u[k])[l]+(D[i,k]+D_reduced[i,k]*u[k] =-grad(X[i,j]+X_reduced[i,j])[j]+Y[i]+Y_reduced[i]*
2520    
2521     M{A} and M{A_reduced} are of rank four, M{B}, M{B_reduced}, M{C} and     *A* and *A_reduced* are of rank four, *B*, *B_reduced*, *C* and
2522     M{C_reduced} are each of rank three, M{D}, M{D_reduced}, M{X_reduced} and     *C_reduced* are each of rank three, *D*, *D_reduced*, *X_reduced* and
2523     M{X} are each of rank two and M{Y} and M{Y_reduced} are of rank one.     *X* are each of rank two and *Y* and *Y_reduced* are of rank one.
2524     The natural boundary conditions take the form:     The natural boundary conditions take the form:
2525    
2526     M{n[j]*((A[i,j,k,l]+A_reduced[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k])+(d[i,k]+d_reduced[i,k])*u[k]=n[j]*(X[i,j]+X_reduced[i,j])+y[i]+y_reduced[i]}     *n[j]*((A[i,j,k,l]+A_reduced[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k])+(d[i,k]+d_reduced[i,k])*u[k]=n[j]*(X[i,j]+X_reduced[i,j])+y[i]+y_reduced[i]*
2527    
2528     The coefficient M{d} is of rank two and M{y} is of rank one both in     The coefficient *d* is of rank two and *y* is of rank one both in
2529     L{FunctionOnBoundary<escript.FunctionOnBoundary>}. The coefficients     `FunctionOnBoundary`. The coefficients
2530     M{d_reduced} is of rank two and M{y_reduced} is of rank one both in     *d_reduced* is of rank two and *y_reduced* is of rank one both in
2531     L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.     `ReducedFunctionOnBoundary`.
2532    
2533     Constraints take the form     Constraints take the form
2534    
2535     M{u[i]=r[i]}  where  M{q[i]>0}     *u[i]=r[i]*  where  *q[i]>0*
2536    
2537     M{r} and M{q} are each rank one. Notice that at some locations not     *r* and *q* are each rank one. Notice that at some locations not
2538     necessarily all components must have a constraint.     necessarily all components must have a constraint.
2539    
2540     The system of PDEs is symmetrical if     The system of PDEs is symmetrical if
2541    
2542        - M{A[i,j,k,l]=A[k,l,i,j]}        - *A[i,j,k,l]=A[k,l,i,j]*
2543        - M{A_reduced[i,j,k,l]=A_reduced[k,l,i,j]}        - *A_reduced[i,j,k,l]=A_reduced[k,l,i,j]*
2544        - M{B[i,j,k]=C[k,i,j]}        - *B[i,j,k]=C[k,i,j]*
2545        - M{B_reduced[i,j,k]=C_reduced[k,i,j]}        - *B_reduced[i,j,k]=C_reduced[k,i,j]*
2546        - M{D[i,k]=D[i,k]}        - *D[i,k]=D[i,k]*
2547        - M{D_reduced[i,k]=D_reduced[i,k]}        - *D_reduced[i,k]=D_reduced[i,k]*
2548        - M{d[i,k]=d[k,i]}        - *d[i,k]=d[k,i]*
2549        - M{d_reduced[i,k]=d_reduced[k,i]}        - *d_reduced[i,k]=d_reduced[k,i]*
2550    
2551     L{LinearPDE} also supports solution discontinuities over a contact region     `LinearPDE` also supports solution discontinuities over a contact region
2552     in the domain. To specify the conditions across the discontinuity we are     in the domain. To specify the conditions across the discontinuity we are
2553     using the generalised flux M{J} which, in the case of a system of PDEs     using the generalised flux *J* which, in the case of a system of PDEs
2554     and several components of the solution, is defined as     and several components of the solution, is defined as
2555    
2556     M{J[i,j]=(A[i,j,k,l]+A_reduced[[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k]-X[i,j]-X_reduced[i,j]}     *J[i,j]=(A[i,j,k,l]+A_reduced[[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k]-X[i,j]-X_reduced[i,j]*
2557    
2558     For the case of single solution component and single PDE M{J} is defined as     For the case of single solution component and single PDE *J* is defined as
2559    
2560     M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[j]+(B[i]+B_reduced[i])*u-X[i]-X_reduced[i]}     *J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[j]+(B[i]+B_reduced[i])*u-X[i]-X_reduced[i]*
2561    
2562     In the context of discontinuities M{n} denotes the normal on the     In the context of discontinuities *n* denotes the normal on the
2563     discontinuity pointing from side 0 towards side 1 calculated from     discontinuity pointing from side 0 towards side 1 calculated from
2564     L{getNormal<escript.FunctionSpace.getNormal>} of L{FunctionOnContactZero<escript.FunctionOnContactZero>}.     `FunctionSpace.getNormal` of `FunctionOnContactZero`.
2565     For a system of PDEs the contact condition takes the form     For a system of PDEs the contact condition takes the form
2566    
2567     M{n[j]*J0[i,j]=n[j]*J1[i,j]=(y_contact[i]+y_contact_reduced[i])- (d_contact[i,k]+d_contact_reduced[i,k])*jump(u)[k]}     *n[j]*J0[i,j]=n[j]*J1[i,j]=(y_contact[i]+y_contact_reduced[i])- (d_contact[i,k]+d_contact_reduced[i,k])*jump(u)[k]*
2568    
2569     where M{J0} and M{J1} are the fluxes on side 0 and side 1 of the     where *J0* and *J1* are the fluxes on side 0 and side 1 of the
2570     discontinuity, respectively. M{jump(u)}, which is the difference of the     discontinuity, respectively. *jump(u)*, which is the difference of the
2571     solution at side 1 and at side 0, denotes the jump of M{u} across     solution at side 1 and at side 0, denotes the jump of *u* across
2572     discontinuity along the normal calculated by L{jump<util.jump>}.     discontinuity along the normal calculated by `jump`.
2573     The coefficient M{d_contact} is of rank two and M{y_contact} is of rank one     The coefficient *d_contact* is of rank two and *y_contact* is of rank one
2574     both in L{FunctionOnContactZero<escript.FunctionOnContactZero>} or     both in `FunctionOnContactZero` or
2575     L{FunctionOnContactOne<escript.FunctionOnContactOne>}.     `FunctionOnContactOne`.
2576     The coefficient M{d_contact_reduced} is of rank two and M{y_contact_reduced}     The coefficient *d_contact_reduced* is of rank two and *y_contact_reduced*
2577     is of rank one both in L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>}     is of rank one both in `ReducedFunctionOnContactZero`
2578     or L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}.     or `ReducedFunctionOnContactOne`.
2579     In case of a single PDE and a single component solution the contact     In case of a single PDE and a single component solution the contact
2580     condition takes the form     condition takes the form
2581    
2582     M{n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)}     *n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)*
2583    
2584     In this case the coefficient M{d_contact} and M{y_contact} are each scalar     In this case the coefficient *d_contact* and *y_contact* are each scalar
2585     both in L{FunctionOnContactZero<escript.FunctionOnContactZero>} or     both in `FunctionOnContactZero` or
2586     L{FunctionOnContactOne<escript.FunctionOnContactOne>} and the coefficient     `FunctionOnContactOne` and the coefficient
2587     M{d_contact_reduced} and M{y_contact_reduced} are each scalar both in     *d_contact_reduced* and *y_contact_reduced* are each scalar both in
2588     L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or     `ReducedFunctionOnContactZero` or
2589     L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}.     `ReducedFunctionOnContactOne`.
2590    
2591     Typical usage::     Typical usage::
2592    
# Line 1628  class LinearPDE(LinearProblem): Line 2600  class LinearPDE(LinearProblem):
2600       """       """
2601       Initializes a new linear PDE.       Initializes a new linear PDE.
2602    
2603       @param domain: domain of the PDE       :param domain: domain of the PDE
2604       @type domain: L{Domain<escript.Domain>}       :type domain: `Domain`
2605       @param numEquations: number of equations. If C{None} the number of       :param numEquations: number of equations. If ``None`` the number of
2606                            equations is extracted from the PDE coefficients.                            equations is extracted from the PDE coefficients.
2607       @param numSolutions: number of solution components. If C{None} the number       :param numSolutions: number of solution components. If ``None`` the number
2608                            of solution components is extracted from the PDE                            of solution components is extracted from the PDE
2609                            coefficients.                            coefficients.
2610       @param debug: if True debug information is printed       :param debug: if True debug information is printed
2611    
2612       """       """
2613       super(LinearPDE, self).__init__(domain,numEquations,numSolutions,debug)       super(LinearPDE, self).__init__(domain,numEquations,numSolutions,debug)
# Line 1670  class LinearPDE(LinearProblem): Line 2642  class LinearPDE(LinearProblem):
2642       """       """
2643       Returns the string representation of the PDE.       Returns the string representation of the PDE.
2644    
2645       @return: a simple representation of the PDE       :return: a simple representation of the PDE
2646       @rtype: C{str}       :rtype: ``str``
2647       """       """
2648       return "<LinearPDE %d>"%id(self)       return "<LinearPDE %d>"%id(self)
2649    
2650     def getRequiredSystemType(self):     def getRequiredOperatorType(self):
2651        """        """
2652        Returns the system type which needs to be used by the current set up.        Returns the system type which needs to be used by the current set up.
2653        """        """
2654        return self.getDomain().getSystemMatrixTypeId(self.getSolverMethod()[0],self.getSolverMethod()[1],self.getSolverPackage(),self.isSymmetric())        if self.isUsingLumping():
2655         return "__ESCRIPT_DATA"
2656          else:
2657         solver_options=self.getSolverOptions()
2658         return self.getDomain().getSystemMatrixTypeId(solver_options.getSolverMethod(), solver_options.getPreconditioner(),solver_options.getPackage(), solver_options.isSymmetric())
2659    
2660     def checkSymmetry(self,verbose=True):     def checkSymmetry(self,verbose=True):
2661        """        """
2662        Tests the PDE for symmetry.        Tests the PDE for symmetry.
2663    
2664        @param verbose: if set to True or not present a report on coefficients        :param verbose: if set to True or not present a report on coefficients
2665                        which break the symmetry is printed.                        which break the symmetry is printed.
2666        @type verbose: C{bool}        :type verbose: ``bool``
2667        @return: True if the PDE is symmetric        :return: True if the PDE is symmetric
2668        @rtype: L{bool}        :rtype: `bool`
2669        @note: This is a very expensive operation. It should be used for        :note: This is a very expensive operation. It should be used for
2670               degugging only! The symmetry flag is not altered.               degugging only! The symmetry flag is not altered.
2671        """        """
2672        out=True        out=True
# Line 1710  class LinearPDE(LinearProblem): Line 2686  class LinearPDE(LinearProblem):
2686         """         """
2687         Returns an instance of a new operator.         Returns an instance of a new operator.
2688         """         """
2689         self.trace("New operator is allocated.")         optype=self.getRequiredOperatorType()
2690           self.trace("New operator of type %s is allocated."%optype)
2691         return self.getDomain().newOperator( \         return self.getDomain().newOperator( \
2692                             self.getNumEquations(), \                             self.getNumEquations(), \
2693                             self.getFunctionSpaceForEquation(), \                             self.getFunctionSpaceForEquation(), \
2694                             self.getNumSolutions(), \                             self.getNumSolutions(), \
2695                             self.getFunctionSpaceForSolution(), \                             self.getFunctionSpaceForSolution(), \
2696                             self.getSystemType())                             optype)
2697    
2698     def getSolution(self,**options):     def getSolution(self):
2699         """         """
2700         Returns the solution of the PDE.         Returns the solution of the PDE.
2701    
2702         @return: the solution         :return: the solution
2703         @rtype: L{Data<escript.Data>}         :rtype: `Data`
        @param options: solver options  
        @keyword verbose: True to get some information during PDE solution  
        @type verbose: C{bool}  
        @keyword reordering: reordering scheme to be used during elimination.  
                             Allowed values are L{NO_REORDERING},  
                             L{MINIMUM_FILL_IN} and L{NESTED_DISSECTION}  
        @keyword iter_max: maximum number of iteration steps allowed  
        @keyword drop_tolerance: threshold for dropping in L{ILUT}  
        @keyword drop_storage: maximum of allowed memory in L{ILUT}  
        @keyword truncation: maximum number of residuals in L{GMRES}  
        @keyword restart: restart cycle length in L{GMRES}  
2704         """         """
2705           option_class=self.getSolverOptions()
2706         if not self.isSolutionValid():         if not self.isSolutionValid():
2707            mat,f=self.getSystem()            mat,f=self.getSystem()
2708            if self.isUsingLumping():            if self.isUsingLumping():
2709             if not util.inf(abs(mat)) > 0.:
2710             raise ZeroDivisionError,"Lumped mass matrix as zero entry (try order 1 elements or HRZ lumping)."
2711               self.setSolution(f*1/mat)               self.setSolution(f*1/mat)
2712            else:            else:
              options[self.TOLERANCE_KEY]=self.getTolerance()  
              options[self.METHOD_KEY]=self.getSolverMethod()[0]  
              options[self.PRECONDITIONER_KEY]=self.getSolverMethod()[1]  
              options[self.PACKAGE_KEY]=self.getSolverPackage()  
              options[self.SYMMETRY_KEY]=self.isSymmetric()  
2713               self.trace("PDE is resolved.")               self.trace("PDE is resolved.")
2714               self.trace("solver options: %s"%str(options))               self.trace("solver options: %s"%str(option_class))
2715               self.setSolution(mat.solve(f,options))               self.setSolution(mat.solve(f,option_class))
2716         return self.getCurrentSolution()         return self.getCurrentSolution()
2717    
2718     def getSystem(self):     def getSystem(self):
2719         """         """
2720         Returns the operator and right hand side of the PDE.         Returns the operator and right hand side of the PDE.
2721    
2722         @return: the discrete version of the PDE         :return: the discrete version of the PDE
2723         @rtype: C{tuple} of L{Operator,<escript.Operator>} and         :rtype: ``tuple`` of `Operator` and
2724                 L{Data<escript.Data>}                 `Data`
2725         """         """
2726         if not self.isOperatorValid() or not self.isRightHandSideValid():         if not self.isOperatorValid() or not self.isRightHandSideValid():
2727            if self.isUsingLumping():            if self.isUsingLumping():
# Line 1786  class LinearPDE(LinearProblem): Line 2750  class LinearPDE(LinearProblem):
2750                   d_reduced=self.getCoefficient("d_reduced")                   d_reduced=self.getCoefficient("d_reduced")
2751                   if not D.isEmpty():                   if not D.isEmpty():
2752                       if self.getNumSolutions()>1:                       if self.getNumSolutions()>1:
2753                          D_times_e=util.matrix_mult(D,numarray.ones((self.getNumSolutions(),)))                          D_times_e=util.matrix_mult(D,numpy.ones((self.getNumSolutions(),)))
2754                       else:                       else:
2755                          D_times_e=D                          D_times_e=D
2756                   else:                   else:
2757                      D_times_e=escript.Data()                      D_times_e=escript.Data()
2758                   if not d.isEmpty():                   if not d.isEmpty():
2759                       if self.getNumSolutions()>1:                       if self.getNumSolutions()>1:
2760                          d_times_e=util.matrix_mult(d,numarray.ones((self.getNumSolutions(),)))                          d_times_e=util.matrix_mult(d,numpy.ones((self.getNumSolutions(),)))
2761                       else:                       else:
2762                          d_times_e=d                          d_times_e=d
2763                   else:                   else:
# Line 1801  class LinearPDE(LinearProblem): Line 2765  class LinearPDE(LinearProblem):
2765    
2766                   if not D_reduced.isEmpty():                   if not D_reduced.isEmpty():
2767                       if self.getNumSolutions()>1:                       if self.getNumSolutions()>1:
2768                          D_reduced_times_e=util.matrix_mult(D_reduced,numarray.ones((self.getNumSolutions(),)))                          D_reduced_times_e=util.matrix_mult(D_reduced,numpy.ones((self.getNumSolutions(),)))
2769                       else:                       else:
2770                          D_reduced_times_e=D_reduced                          D_reduced_times_e=D_reduced
2771                   else:                   else:
2772                      D_reduced_times_e=escript.Data()                      D_reduced_times_e=escript.Data()
2773                   if not d_reduced.isEmpty():                   if not d_reduced.isEmpty():
2774                       if self.getNumSolutions()>1:                       if self.getNumSolutions()>1:
2775                          d_reduced_times_e=util.matrix_mult(d_reduced,numarray.ones((self.getNumSolutions(),)))                          d_reduced_times_e=util.matrix_mult(d_reduced,numpy.ones((self.getNumSolutions(),)))
2776                       else:                       else:
2777                          d_reduced_times_e=d_reduced                          d_reduced_times_e=d_reduced
2778                   else:                   else:
# Line 1816  class LinearPDE(LinearProblem): Line 2780  class LinearPDE(LinearProblem):
2780    
2781                   self.resetOperator()                   self.resetOperator()
2782                   operator=self.getCurrentOperator()                   operator=self.getCurrentOperator()
2783                   if False and hasattr(self.getDomain(), "addPDEToLumpedSystem") :                   if hasattr(self.getDomain(), "addPDEToLumpedSystem") :
2784                      self.getDomain().addPDEToLumpedSystem(operator, D_times_e, d_times_e)              hrz_lumping=( self.getSolverOptions().getSolverMethod() ==  SolverOptions.HRZ_LUMPING )
2785                      self.getDomain().addPDEToLumpedSystem(operator, D_reduced_times_e, d_reduced_times_e)              self.getDomain().addPDEToLumpedSystem(operator, D_times_e, d_times_e,  hrz_lumping )
2786                self.getDomain().addPDEToLumpedSystem(operator, D_reduced_times_e, d_reduced_times_e, hrz_lumping)
2787                   else:                   else:
2788                      self.getDomain().addPDEToRHS(operator, \                      self.getDomain().addPDEToRHS(operator, \
2789                                                   escript.Data(), \                                                   escript.Data(), \
# Line 1924  class LinearPDE(LinearProblem): Line 2889  class LinearPDE(LinearProblem):
2889                   self.insertConstraint(rhs_only=False)                   self.insertConstraint(rhs_only=False)
2890                   self.trace("New operator has been built.")                   self.trace("New operator has been built.")
2891                   self.validOperator()                   self.validOperator()
2892           self.setSystemStatus()
2893           self.trace("System status is %s."%self.getSystemStatus())
2894         return (self.getCurrentOperator(), self.getCurrentRightHandSide())         return (self.getCurrentOperator(), self.getCurrentRightHandSide())
2895    
2896     def insertConstraint(self, rhs_only=False):     def insertConstraint(self, rhs_only=False):
2897        """        """
2898        Applies the constraints defined by q and r to the PDE.        Applies the constraints defined by q and r to the PDE.
2899    
2900        @param rhs_only: if True only the right hand side is altered by the        :param rhs_only: if True only the right hand side is altered by the
2901                         constraint                         constraint
2902        @type rhs_only: C{bool}        :type rhs_only: ``bool``
2903        """        """
2904        q=self.getCoefficient("q")        q=self.getCoefficient("q")
2905        r=self.getCoefficient("r")        r=self.getCoefficient("r")
# Line 1960  class LinearPDE(LinearProblem): Line 2927  class LinearPDE(LinearProblem):
2927        """        """
2928        Sets new values to coefficients.        Sets new values to coefficients.
2929    
2930        @param coefficients: new values assigned to coefficients        :param coefficients: new values assigned to coefficients
2931        @keyword A: value for coefficient C{A}        :keyword A: value for coefficient ``A``
2932        @type A: any type that can be cast to a L{Data<escript.Data>} object on        :type A: any type that can be cast to a `Data` object on
2933                 L{Function<escript.Function>}                 `Function`
2934        @keyword A_reduced: value for coefficient C{A_reduced}        :keyword A_reduced: value for coefficient ``A_reduced``
2935        @type A_reduced: any type that can be cast to a L{Data<escript.Data>}        :type A_reduced: any type that can be cast to a `Data`
2936                         object on L{ReducedFunction<escript.ReducedFunction>}                         object on `ReducedFunction`
2937        @keyword B: value for coefficient C{B}        :keyword B: value for coefficient ``B``
2938        @type B: any type that can be cast to a L{Data<escript.Data>} object on        :type B: any type that can be cast to a `Data` object on
2939                 L{Function<escript.Function>}                 `Function`
2940        @keyword B_reduced: value for coefficient C{B_reduced}        :keyword B_reduced: value for coefficient ``B_reduced``
2941        @type B_reduced: any type that can be cast to a L{Data<escript.Data>}        :type B_reduced: any type that can be cast to a `Data`
2942                         object on L{ReducedFunction<escript.ReducedFunction>}                         object on `ReducedFunction`
2943        @keyword C: value for coefficient C{C}        :keyword C: value for coefficient ``C``
2944        @type C: any type that can be cast to a L{Data<escript.Data>} object on        :type C: any type that can be cast to a `Data` object on
2945                 L{Function<escript.Function>}                 `Function`
2946        @keyword C_reduced: value for coefficient C{C_reduced}        :keyword C_reduced: value for coefficient ``C_reduced``
2947        @type C_reduced: any type that can be cast to a L{Data<escript.Data>}        :type C_reduced: any type that can be cast to a `Data`
2948                         object on L{ReducedFunction<escript.ReducedFunction>}                         object on `ReducedFunction`
2949        @keyword D: value for coefficient C{D}        :keyword D: value for coefficient ``D``
2950        @type D: any type that can be cast to a L{Data<escript.Data>} object on        :type D: any type that can be cast to a `Data` object on
2951                 L{Function<escript.Function>}                 `Function`
2952        @keyword D_reduced: value for coefficient C{D_reduced}        :keyword D_reduced: value for coefficient ``D_reduced``
2953        @type D_reduced: any type that can be cast to a L{Data<escript.Data>}        :type D_reduced: any type that can be cast to a `Data`
2954                         object on L{ReducedFunction<escript.ReducedFunction>}                         object on `ReducedFunction`
2955        @keyword X: value for coefficient C{X}        :keyword X: value for coefficient ``X``
2956        @type X: any type that can be cast to a L{Data<escript.Data>} object on        :type X: any type that can be cast to a `Data` object on
2957                 L{Function<escript.Function>}                 `Function`
2958        @keyword X_reduced: value for coefficient C{X_reduced}        :keyword X_reduced: value for coefficient ``X_reduced``
2959        @type X_reduced: any type that can be cast to a L{Data<escript.Data>}        :type X_reduced: any type that can be cast to a `Data`
2960                         object on L{ReducedFunction<escript.ReducedFunction>}                         object on `ReducedFunction`
2961        @keyword Y: value for coefficient C{Y}        :keyword Y: value for coefficient ``Y``
2962        @type Y: any type that can be cast to a L{Data<escript.Data>} object on        :type Y: any type that can be cast to a `Data` object on
2963                 L{Function<escript.Function>}                 `Function`
2964        @keyword Y_reduced: value for coefficient C{Y_reduced}        :keyword Y_reduced: value for coefficient ``Y_reduced``
2965        @type Y_reduced: any type that can be cast to a L{Data<escript.Data>}        :type Y_reduced: any type that can be cast to a `Data`
2966                         object on L{ReducedFunction<escript.Function>}                         object on `ReducedFunction`
2967        @keyword d: value for coefficient C{d}        :keyword d: value for coefficient ``d``
2968        @type d: any type that can be cast to a L{Data<escript.Data>} object on        :type d: any type that can be cast to a `Data` object on
2969                 L{FunctionOnBoundary<escript.FunctionOnBoundary>}                 `FunctionOnBoundary`
2970        @keyword d_reduced: value for coefficient C{d_reduced}        :keyword d_reduced: value for coefficient ``d_reduced``
2971        @type d_reduced: any type that can be cast to a L{Data<escript.Data>}        :type d_reduced: any type that can be cast to a `Data`
2972                         object on L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}                         object on `ReducedFunctionOnBoundary`
2973        @keyword y: value for coefficient C{y}        :keyword y: value for coefficient ``y``
2974        @type y: any type that can be cast to a L{Data<escript.Data>} object on        :type y: any type that can be cast to a `Data` object on
2975                 L{FunctionOnBoundary<escript.FunctionOnBoundary>}                 `FunctionOnBoundary`
2976        @keyword d_contact: value for coefficient C{d_contact}        :keyword d_contact: value for coefficient ``d_contact``
2977        @type d_contact: any type that can be cast to a L{Data<escript.Data>}        :type d_contact: any type that can be cast to a `Data`
2978                         object on L{FunctionOnContactOne<escript.FunctionOnContactOne>}                         object on `FunctionOnContactOne`
2979                         or L{FunctionOnContactZero<escript.FunctionOnContactZero>}                         or `FunctionOnContactZero`
2980        @keyword d_contact_reduced: value for coefficient C{d_contact_reduced}        :keyword d_contact_reduced: value for coefficient ``d_contact_reduced``
2981        @type d_contact_reduced: any type that can be cast to a L{Data<escript.Data>}        :type d_contact_reduced: any type that can be cast to a `Data`
2982                                 object on L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}                                 object on `ReducedFunctionOnContactOne`
2983                                 or L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>}                                 or `ReducedFunctionOnContactZero`
2984        @keyword y_contact: value for coefficient C{y_contact}        :keyword y_contact: value for coefficient ``y_contact``
2985        @type y_contact: any type that can be cast to a L{Data<escript.Data>}        :type y_contact: any type that can be cast to a `Data`
2986                         object on L{FunctionOnContactOne<escript.FunctionOnContactOne>}                         object on `FunctionOnContactOne`
2987                         or L{FunctionOnContactZero<escript.FunctionOnContactZero>}                         or `FunctionOnContactZero`
2988        @keyword y_contact_reduced: value for coefficient C{y_contact_reduced}        :keyword y_contact_reduced: value for coefficient ``y_contact_reduced``
2989        @type y_contact_reduced: any type that can be cast to a L{Data<escript.Data>}        :type y_contact_reduced: any type that can be cast to a `Data`
2990                                 object on L{ReducedFunctionOnContactOne<escript.FunctionOnContactOne>}                                 object on `ReducedFunctionOnContactOne`
2991                                 or L{ReducedFunctionOnContactZero<escript.FunctionOnContactZero>}                                 or `ReducedFunctionOnContactZero`
2992        @keyword r: values prescribed to the solution at the locations of        :keyword r: values prescribed to the solution at the locations of
2993                    constraints                    constraints
2994        @type r: any type that can be cast to a L{Data<escript.Data>} object on        :type r: any type that can be cast to a `Data` object on
2995                 L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}                 `Solution` or `ReducedSolution`
2996                 depending on whether reduced order is used for the solution                 depending on whether reduced order is used for the solution
2997        @keyword q: mask for location of constraints        :keyword q: mask for location of constraints
2998        @type q: any type that can be cast to a L{Data<escript.Data>} object on        :type q: any type that can be cast to a `Data` object on
2999                 L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}                 `Solution` or `ReducedSolution`
3000                 depending on whether reduced order is used for the                 depending on whether reduced order is used for the
3001                 representation of the equation                 representation of the equation
3002        @raise IllegalCoefficient: if an unknown coefficient keyword is used        :raise IllegalCoefficient: if an unknown coefficient keyword is used
3003        """        """
3004        super(LinearPDE,self).setValue(**coefficients)        super(LinearPDE,self).setValue(**coefficients)
3005        # check if the systrem is inhomogeneous:        # check if the systrem is inhomogeneous:
# Line 2049  class LinearPDE(LinearProblem): Line 3016  class LinearPDE(LinearProblem):
3016       """       """
3017       Returns the residual of u or the current solution if u is not present.       Returns the residual of u or the current solution if u is not present.
3018    
3019       @param u: argument in the residual calculation. It must be representable       :param u: argument in the residual calculation. It must be representable
3020                 in L{self.getFunctionSpaceForSolution()}. If u is not present                 in `self.getFunctionSpaceForSolution()`. If u is not present
3021                 or equals C{None} the current solution is used.                 or equals ``None`` the current solution is used.
3022       @type u: L{Data<escript.Data>} or None       :type u: `Data` or None
3023       @return: residual of u       :return: residual of u
3024       @rtype: L{Data<escript.Data>}       :rtype: `Data`
3025       """       """
3026       if u==None:       if u==None:
3027          return self.getOperator()*self.getSolution()-self.getRightHandSide()          return self.getOperator()*self.getSolution()-self.getRightHandSide()
# Line 2063  class LinearPDE(LinearProblem): Line 3030  class LinearPDE(LinearProblem):
3030    
3031     def getFlux(self,u=None):     def getFlux(self,u=None):
3032       """       """
3033       Returns the flux M{J} for a given M{u}.       Returns the flux *J* for a given *u*.
3034    
3035       M{J[i,j]=(A[i,j,k,l]+A_reduced[A[i,j,k,l]]*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])u[k]-X[i,j]-X_reduced[i,j]}       *J[i,j]=(A[i,j,k,l]+A_reduced[A[i,j,k,l]]*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])u[k]-X[i,j]-X_reduced[i,j]*
3036    
3037       or       or
3038    
3039       M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[l]+(B[j]+B_reduced[j])u-X[j]-X_reduced[j]}       *J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[l]+(B[j]+B_reduced[j])u-X[j]-X_reduced[j]*
3040    
3041       @param u: argument in the flux. If u is not present or equals L{None} the       :param u: argument in the flux. If u is not present or equals `None` the
3042                 current solution is used.                 current solution is used.
3043       @type u: L{Data<escript.Data>} or None       :type u: `Data` or None
3044       @return: flux       :return: flux
3045       @rtype: L{Data<escript.Data>}       :rtype: `Data`
3046       """       """
3047       if u==None: u=self.getSolution()       if u==None: u=self.getSolution()
3048       return util.tensormult(self.getCoefficient("A"),util.grad(u,Funtion(self.getDomain))) \       return util.tensormult(self.getCoefficient("A"),util.grad(u,Funtion(self.getDomain))) \
# Line 2089  class LinearPDE(LinearProblem): Line 3056  class LinearPDE(LinearProblem):
3056  class Poisson(LinearPDE):  class Poisson(LinearPDE):
3057     """     """
3058     Class to define a Poisson equation problem. This is generally a     Class to define a Poisson equation problem. This is generally a
3059     L{LinearPDE} of the form     `LinearPDE` of the form
3060    
3061     M{-grad(grad(u)[j])[j] = f}     *-grad(grad(u)[j])[j] = f*
3062    
3063     with natural boundary conditions     with natural boundary conditions
3064    
3065     M{n[j]*grad(u)[j] = 0 }     *n[j]*grad(u)[j] = 0*
3066    
3067     and constraints:     and constraints:
3068    
3069     M{u=0} where M{q>0}     *u=0* where *q>0*
3070    
3071     """     """
3072    
# Line 2107  class Poisson(LinearPDE): Line 3074  class Poisson(LinearPDE):
3074       """       """
3075       Initializes a new Poisson equation.       Initializes a new Poisson equation.
3076    
3077       @param domain: domain of the PDE       :param domain: domain of the PDE
3078       @type domain: L{Domain<escript.Domain>}       :type domain: `Domain`
3079       @param debug: if True debug information is printed       :param debug: if True debug information is printed
3080    
3081       """       """
3082       super(Poisson, self).__init__(domain,1,1,debug)       super(Poisson, self).__init__(domain,1,1,debug)
# Line 2122  class Poisson(LinearPDE): Line 3089  class Poisson(LinearPDE):
3089       """       """
3090       Sets new values to coefficients.       Sets new values to coefficients.
3091    
3092       @param coefficients: new values assigned to coefficients       :param coefficients: new values assigned to coefficients
3093       @keyword f: value for right hand side M{f}       :keyword f: value for right hand side *f*
3094       @type f: any type that can be cast to a L{Scalar<escript.Scalar>} object       :type f: any type that can be cast to a `Scalar` object
3095                on L{Function<escript.Function>}                on `Function`
3096       @keyword q: mask for location of constraints       :keyword q: mask for location of constraints
3097       @type q: any type that can be cast to a rank zero L{Data<escript.Data>}       :type q: any type that can be cast to a rank zero `Data`
3098                object on L{Solution<escript.Solution>} or                object on `Solution` or
3099                L{ReducedSolution<escript.ReducedSolution>} depending on whether                `ReducedSolution` depending on whether
3100                reduced order is used for the representation of the equation                reduced order is used for the representation of the equation
3101       @raise IllegalCoefficient: if an unknown coefficient keyword is used       :raise IllegalCoefficient: if an unknown coefficient keyword is used
3102       """       """
3103       super(Poisson, self).setValue(**coefficients)       super(Poisson, self).setValue(**coefficients)
3104    
3105    
3106     def getCoefficient(self,name):     def getCoefficient(self,name):
3107       """       """
3108       Returns the value of the coefficient C{name} of the general PDE.       Returns the value of the coefficient ``name`` of the general PDE.
3109    
3110       @param name: name of the coefficient requested       :param name: name of the coefficient requested
3111       @type name: C{string}       :type name: ``string``
3112       @return: the value of the coefficient C{name}       :return: the value of the coefficient ``name``
3113       @rtype: L{Data<escript.Data>}       :rtype: `Data`
3114       @raise IllegalCoefficient: invalid coefficient name       :raise IllegalCoefficient: invalid coefficient name
3115       @note: This method is called by the assembling routine to map the Poisson       :note: This method is called by the assembling routine to map the Poisson
3116              equation onto the general PDE.              equation onto the general PDE.
3117       """       """
3118       if name == "A" :       if name == "A" :
# Line 2160  class Poisson(LinearPDE): Line 3127  class Poisson(LinearPDE):
3127  class Helmholtz(LinearPDE):  class Helmholtz(LinearPDE):
3128     """     """
3129     Class to define a Helmholtz equation problem. This is generally a     Class to define a Helmholtz equation problem. This is generally a
3130     L{LinearPDE} of the form     `LinearPDE` of the form
3131    
3132     M{S{omega}*u - grad(k*grad(u)[j])[j] = f}     *omega*u - grad(k*grad(u)[j])[j] = f*
3133    
3134     with natural boundary conditions     with natural boundary conditions
3135    
3136     M{k*n[j]*grad(u)[j] = g- S{alpha}u }     *k*n[j]*grad(u)[j] = g- alphau*
3137    
3138     and constraints:     and constraints:
3139    
3140     M{u=r} where M{q>0}     *u=r* where *q>0*
3141    
3142     """     """
3143    
# Line 2178  class Helmholtz(LinearPDE): Line 3145  class Helmholtz(LinearPDE):
3145       """       """
3146       Initializes a new Helmholtz equation.       Initializes a new Helmholtz equation.
3147    
3148       @param domain: domain of the PDE       :param domain: domain of the PDE
3149       @type domain: L{Domain<escript.Domain>}       :type domain: `Domain`
3150       @param debug: if True debug information is printed       :param debug: if True debug information is printed
3151    
3152       """       """
3153       super(Helmholtz, self).__init__(domain,1,1,debug)       super(Helmholtz, self).__init__(domain,1,1,debug)
# Line 2198  class Helmholtz(LinearPDE): Line 3165  class Helmholtz(LinearPDE):
3165       """       """
3166       Sets new values to coefficients.       Sets new values to coefficients.
3167    
3168       @param coefficients: new values assigned to coefficients       :param coefficients: new values assigned to coefficients
3169       @keyword omega: value for coefficient M{S{omega}}       :keyword omega: value for coefficient *omega*
3170       @type omega: any type that can be cast to a L{Scalar<escript.Scalar>}       :type omega: any type that can be cast to a `Scalar`
3171                    object on L{Function<escript.Function>}                    object on `Function`
3172       @keyword k: value for coefficient M{k}       :keyword k: value for coefficient *k*
3173       @type k: any type that can be cast to a L{Scalar<escript.Scalar>} object       :type k: any type that can be cast to a `Scalar` object
3174                on L{Function<escript.Function>}                on `Function`
3175       @keyword f: value for right hand side M{f}       :keyword f: value for right hand side *f*
3176       @type f: any type that can be cast to a L{Scalar<escript.Scalar>} object       :type f: any type that can be cast to a `Scalar` object
3177                on L{Function<escript.Function>}                on `Function`
3178       @keyword alpha: value for right hand side M{S{alpha}}       :keyword alpha: value for right hand side *alpha*
3179       @type alpha: any type that can be cast to a L{Scalar<escript.Scalar>}       :type alpha: any type that can be cast to a `Scalar`
3180                    object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}                    object on `FunctionOnBoundary`
3181       @keyword g: value for right hand side M{g}       :keyword g: value for right hand side *g*
3182       @type g: any type that can be cast to a L{Scalar<escript.Scalar>} object       :type g: any type that can be cast to a `Scalar` object
3183                on L{FunctionOnBoundary<escript.FunctionOnBoundary>}                on `FunctionOnBoundary`
3184       @keyword r: prescribed values M{r} for the solution in constraints       :keyword r: prescribed values *r* for the solution in constraints
3185       @type r: any type that can be cast to a L{Scalar<escript.Scalar>} object       :type r: any type that can be cast to a `Scalar` object
3186                on L{Solution<escript.Solution>} or                on `Solution` or
3187                L{ReducedSolution<escript.ReducedSolution>} depending on whether                `ReducedSolution` depending on whether
3188                reduced order is used for the representation of the equation                reduced order is used for the representation of the equation
3189       @keyword q: mask for the location of constraints       :keyword q: mask for the location of constraints
3190       @type q: any type that can be cast to a L{Scalar<escript.Scalar>} object       :type q: any type that can be cast to a `Scalar` object
3191                on L{Solution<escript.Solution>} or                on `Solution` or
3192                L{ReducedSolution<escript.ReducedSolution>} depending on whether                `ReducedSolution` depending on whether
3193                reduced order is used for the representation of the equation                reduced order is used for the representation of the equation
3194       @raise IllegalCoefficient: if an unknown coefficient keyword is used       :raise IllegalCoefficient: if an unknown coefficient keyword is used
3195       """       """
3196       super(Helmholtz, self).setValue(**coefficients)       super(Helmholtz, self).setValue(**coefficients)
3197    
3198     def getCoefficient(self,name):     def getCoefficient(self,name):
3199       """       """
3200       Returns the value of the coefficient C{name} of the general PDE.       Returns the value of the coefficient ``name`` of the general PDE.
3201    
3202       @param name: name of the coefficient requested       :param name: name of the coefficient requested
3203       @type name: C{string}       :type name: ``string``
3204       @return: the value of the coefficient C{name}       :return: the value of the coefficient ``name``
3205       @rtype: L{Data<escript.Data>}       :rtype: `Data`
3206       @raise IllegalCoefficient: invalid name       :raise IllegalCoefficient: invalid name
3207       """       """
3208       if name == "A" :       if name == "A" :
3209           if self.getCoefficient("k").isEmpty():           if self.getCoefficient("k").isEmpty():
3210                return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))                return escript.Data(numpy.identity(self.getDim()),escript.Function(self.getDomain()))
3211           else:           else:
3212                return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))*self.getCoefficient("k")                return escript.Data(numpy.identity(self.getDim()),escript.Function(self.getDomain()))*self.getCoefficient("k")
3213       elif name == "D" :       elif name == "D" :
3214           return self.getCoefficient("omega")           return self.getCoefficient("omega")
3215       elif name == "Y" :       elif name == "Y" :
# Line 2262  class LameEquation(LinearPDE): Line 3229  class LameEquation(LinearPDE):
3229     """     """
3230     Class to define a Lame equation problem. This problem is defined as:     Class to define a Lame equation problem. This problem is defined as:
3231    
3232     M{-grad(S{mu}*(grad(u[i])[j]+grad(u[j])[i]))[j] - grad(S{lambda}*grad(u[k])[k])[j] = F_i -grad(S{sigma}[ij])[j] }     *-grad(mu*(grad(u[i])[j]+grad(u[j])[i]))[j] - grad(lambda*grad(u[k])[k])[j] = F_i -grad(sigma[ij])[j]*
3233    
3234     with natural boundary conditions:     with natural boundary conditions:
3235    
3236     M{n[j]*(S{mu}*(grad(u[i])[j]+grad(u[j])[i]) + S{lambda}*grad(u[k])[k]) = f_i +n[j]*S{sigma}[ij] }     *n[j]*(mu*(grad(u[i])[j]+grad(u[j])[i]) + lambda*grad(u[k])[k]) = f_i +n[j]*sigma[ij]*
3237    
3238     and constraints:     and constraints:
3239    
3240     M{u[i]=r[i]} where M{q[i]>0}     *u[i]=r[i]* where *q[i]>0*
3241    
3242     """     """
3243    
# Line 2278  class LameEquation(LinearPDE): Line 3245  class LameEquation(LinearPDE):
3245        """        """
3246        Initializes a new Lame equation.        Initializes a new Lame equation.
3247    
3248        @param domain: domain of the PDE        :param domain: domain of the PDE
3249        @type domain: L{Domain<escript.Domain>}        :type domain: `Domain`
3250        @param debug: if True debug information is printed        :param debug: if True debug information is printed
3251    
3252        """        """
3253        super(LameEquation, self).__init__(domain,\        super(LameEquation, self).__init__(domain,\
# Line 2296  class LameEquation(LinearPDE): Line 3263  class LameEquation(LinearPDE):
3263       """       """
3264       Sets new values to coefficients.       Sets new values to coefficients.
3265    
3266       @param coefficients: new values assigned to coefficients       :param coefficients: new values assigned to coefficients
3267       @keyword lame_mu: value for coefficient M{S{mu}}       :keyword lame_mu: value for coefficient *mu*
3268       @type lame_mu: any type that can be cast to a L{Scalar<escript.Scalar>}       :type lame_mu: any type that can be cast to a `Scalar`
3269                      object on L{Function<escript.Function>}                      object on `Function`
3270       @keyword lame_lambda: value for coefficient M{S{lambda}}       :keyword lame_lambda: value for coefficient *lambda*
3271       @type lame_lambda: any type that can be cast to a L{Scalar<escript.Scalar>}       :type lame_lambda: any type that can be cast to a `Scalar`
3272                          object on L{Function<escript.Function>}                          object on `Function`
3273       @keyword F: value for internal force M{F}       :keyword F: value for internal force *F*
3274       @type F: any type that can be cast to a L{Vector<escript.Vector>} object       :type F: any type that can be cast to a `Vector` object
3275                on L{Function<escript.Function>}                on `Function`
3276       @keyword sigma: value for initial stress M{S{sigma}}       :keyword sigma: value for initial stress *sigma*
3277       @type sigma: any type that can be cast to a L{Tensor<escript.Tensor>}       :type sigma: any type that can be cast to a `Tensor`
3278                    object on L{Function<escript.Function>}                    object on `Function`
3279       @keyword f: value for external force M{f}       :keyword f: value for external force *f*
3280       @type f: any type that can be cast to a L{Vector<escript.Vector>} object       :type f: any type that can be cast to a `Vector` object
3281                on L{FunctionOnBoundary<escript.FunctionOnBoundary>}                on `FunctionOnBoundary`
3282       @keyword r: prescribed values M{r} for the solution in constraints       :keyword r: prescribed values *r* for the solution in constraints
3283       @type r: any type that can be cast to a L{Vector<escript.Vector>} object       :type r: any type that can be cast to a `Vector` object
3284                on L{Solution<escript.Solution>} or                on `Solution` or
3285                L{ReducedSolution<escript.ReducedSolution>} depending on whether                `ReducedSolution` depending on whether
3286                reduced order is used for the representation of the equation                reduced order is used for the representation of the equation
3287       @keyword q: mask for the location of constraints       :keyword q: mask for the location of constraints
3288       @type q: any type that can be cast to a L{Vector<escript.Vector>} object       :type q: any type that can be cast to a `Vector` object
3289                on L{Solution<escript.Solution>} or                on `Solution` or
3290                L{ReducedSolution<escript.ReducedSolution>} depending on whether                `ReducedSolution` depending on whether
3291                reduced order is used for the representation of the equation                reduced order is used for the representation of the equation
3292       @raise IllegalCoefficient: if an unknown coefficient keyword is used       :raise IllegalCoefficient: if an unknown coefficient keyword is used
3293       """       """
3294       super(LameEquation, self).setValues(**coefficients)       super(LameEquation, self).setValues(**coefficients)
3295    
3296     def getCoefficient(self,name):     def getCoefficient(self,name):
3297       """       """
3298       Returns the value of the coefficient C{name} of the general PDE.       Returns the value of the coefficient ``name`` of the general PDE.
3299    
3300       @param name: name of the coefficient requested       :param name: name of the coefficient requested
3301       @type name: C{string}       :type name: ``string``
3302       @return: the value of the coefficient C{name}       :return: the value of the coefficient ``name``
3303       @rtype: L{Data<escript.Data>}       :rtype: `Data`
3304       @raise IllegalCoefficient: invalid coefficient name       :raise IllegalCoefficient: invalid coefficient name
3305       """       """
3306       out =self.createCoefficient("A")       out =self.createCoefficient("A")
3307       if name == "A" :       if name == "A" :
# Line 2367  class LameEquation(LinearPDE): Line 3334  class LameEquation(LinearPDE):
3334       else:       else:
3335          return super(LameEquation, self).getCoefficient(name)          return super(LameEquation, self).getCoefficient(name)
3336    
3337    
3338  def LinearSinglePDE(domain,debug=False):  def LinearSinglePDE(domain,debug=False):
3339     """     """
3340     Defines a single linear PDE.     Defines a single linear PDE.
3341    
3342     @param domain: domain of the PDE     :param domain: domain of the PDE
3343     @type domain: L{Domain<escript.Domain>}     :type domain: `Domain`
3344     @param debug: if True debug information is printed     :param debug: if True debug information is printed
3345     @rtype: L{LinearPDE}     :rtype: `LinearPDE`
3346     """     """
3347     return LinearPDE(domain,numEquations=1,numSolutions=1,debug=debug)     return LinearPDE(domain,numEquations=1,numSolutions=1,debug=debug)
3348    
# Line 2382  def LinearPDESystem(domain,debug=False): Line 3350  def LinearPDESystem(domain,debug=False):
3350     """     """
3351     Defines a system of linear PDEs.     Defines a system of linear PDEs.
3352    
3353     @param domain: domain of the PDEs     :param domain: domain of the PDEs
3354     @type domain: L{Domain<escript.Domain>}     :type domain: `Domain`
3355     @param debug: if True debug information is printed     :param debug: if True debug information is printed
3356     @rtype: L{LinearPDE}     :rtype: `LinearPDE`
3357     """     """
3358     return LinearPDE(domain,numEquations=domain.getDim(),numSolutions=domain.getDim(),debug=debug)     return LinearPDE(domain,numEquations=domain.getDim(),numSolutions=domain.getDim(),debug=debug)
3359    
# Line 2394  class TransportPDE(LinearProblem): Line 3362  class TransportPDE(LinearProblem):
3362     """     """
3363     This class is used to define a transport problem given by a general linear,     This class is used to define a transport problem given by a general linear,
3364     time dependent, second order PDE for an unknown, non-negative,     time dependent, second order PDE for an unknown, non-negative,
3365     time-dependent function M{u} on a given domain defined through a     time-dependent function *u* on a given domain defined through a
3366     L{Domain<escript.Domain>} object.     `Domain` object.
3367    
3368     For a single equation with a solution with a single component the transport     For a single equation with a solution with a single component the transport
3369     problem is defined in the following form:     problem is defined in the following form:
3370    
3371     M{(M+M_reduced)*u_t=-(grad(A[j,l]+A_reduced[j,l])*grad(u)[l]+(B[j]+B_reduced[j])u)[j]+(C[l]+C_reduced[l])*grad(u)[l]+(D+D_reduced)-grad(X+X_reduced)[j,j]+(Y+Y_reduced)}     *(M+M_reduced)*u_t=-(grad(A[j,l]+A_reduced[j,l]) * grad(u)[l]+(B[j]+B_reduced[j])u)[j]+(C[l]+C_reduced[l])*grad(u)[l]+(D+D_reduced)-grad(X+X_reduced)[j,j]+(Y+Y_reduced)*
3372    
3373     where M{u_t} denotes the time derivative of M{u} and M{grad(F)} denotes the     where *u_t* denotes the time derivative of *u* and *grad(F)* denotes the
3374     spatial derivative of M{F}. Einstein's summation convention,  ie. summation     spatial derivative of *F*. Einstein's summation convention,  ie. summation
3375     over indexes appearing twice in a term of a sum performed, is used.     over indexes appearing twice in a term of a sum performed, is used.
3376     The coefficients M{M}, M{A}, M{B}, M{C}, M{D}, M{X} and M{Y} have to be     The coefficients *M*, *A*, *B*, *C*, *D*, *X* and *Y* have to be
3377     specified through L{Data<escript.Data>} objects in L{Function<escript.Function>}     specified through `Data` objects in `Function`
3378     and the coefficients M{M_reduced}, M{A_reduced}, M{B_reduced}, M{C_reduced},     and the coefficients *M_reduced*, *A_reduced*, *B_reduced*, *C_reduced*,
3379     M{D_reduced}, M{X_reduced} and M{Y_reduced} have to be specified through     *D_reduced*, *X_reduced* and *Y_reduced* have to be specified through
3380     L{Data<escript.Data>} objects in L{ReducedFunction<escript.ReducedFunction>}.     `Data` objects in `ReducedFunction`.
3381     It is also allowed to use objects that can be converted into such     It is also allowed to use objects that can be converted into such
3382     L{Data<escript.Data>} objects. M{M} and M{M_reduced} are scalar, M{A} and     `Data` objects. *M* and *M_reduced* are scalar, *A* and
3383     M{A_reduced} are rank two, M{B}, M{C}, M{X}, M{B_reduced}, M{C_reduced} and     *A_reduced* are rank two, *B*, *C*, *X*, *B_reduced*, *C_reduced* and
3384     M{X_reduced} are rank one and M{D}, M{D_reduced}, M{Y} and M{Y_reduced} are     *X_reduced* are rank one and *D*, *D_reduced*, *Y* and *Y_reduced* are
3385     scalar.     scalar.
3386    
3387     The following natural boundary conditions are considered:     The following natural boundary conditions are considered:
3388    
3389     M{n[j]*((A[i,j]+A_reduced[i,j])*grad(u)[l]+(B+B_reduced)[j]*u+X[j]+X_reduced[j])+(d+d_reduced)*u+y+y_reduced=(m+m_reduced)*u_t}     *n[j]*((A[i,j]+A_reduced[i,j])*grad(u)[l]+(B+B_reduced)[j]*u+X[j]+X_reduced[j])+(d+d_reduced)*u+y+y_reduced=(m+m_reduced)*u_t*
3390    
3391     where M{n} is the outer normal field. Notice that the coefficients M{A},     where *n* is the outer normal field. Notice that the coefficients *A*,
3392     M{A_reduced}, M{B}, M{B_reduced}, M{X} and M{X_reduced} are defined in the     *A_reduced*, *B*, *B_reduced*, *X* and *X_reduced* are defined in the
3393     transport problem. The coefficients M{m}, M{d} and M{y} are each a scalar in     transport problem. The coefficients *m*, *d* and *y* are each a scalar in
3394     L{FunctionOnBoundary<escript.FunctionOnBoundary>} and the coefficients     `FunctionOnBoundary` and the coefficients
3395     M{m_reduced}, M{d_reduced} and M{y_reduced} are each a scalar in     *m_reduced*, *d_reduced* and *y_reduced* are each a scalar in
3396     L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.     `ReducedFunctionOnBoundary`.
3397    
3398     Constraints for the solution prescribing the value of the solution at     Constraints for the solution prescribing the value of the solution at
3399     certain locations in the domain have the form     certain locations in the domain have the form
3400    
3401     M{u_t=r} where M{q>0}     *u_t=r* where *q>0*
3402    
3403     M{r} and M{q} are each scalar where M{q} is the characteristic function     *r* and *q* are each scalar where *q* is the characteristic function
3404     defining where the constraint is applied. The constraints override any other     defining where the constraint is applied. The constraints override any other
3405     condition set by the transport problem or the boundary condition.     condition set by the transport problem or the boundary condition.
3406    
3407     The transport problem is symmetrical if     The transport problem is symmetrical if
3408    
3409     M{A[i,j]=A[j,i]} and M{B[j]=C[j]} and M{A_reduced[i,j]=A_reduced[j,i]} and     *A[i,j]=A[j,i]* and *B[j]=C[j]* and *A_reduced[i,j]=A_reduced[j,i]* and
3410     M{B_reduced[j]=C_reduced[j]}     *B_reduced[j]=C_reduced[j]*
3411    
3412     For a system and a solution with several components the transport problem     For a system and a solution with several components the transport problem
3413     has the form     has the form
3414    
3415     M{(M[i,k]+M_reduced[i,k])*u[k]_t=-grad((A[i,j,k,l]+A_reduced[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k])[j]+(C[i,k,l]+C_reduced[i,k,l])*grad(u[k])[l]+(D[i,k]+D_reduced[i,k]*u[k]-grad(X[i,j]+X_reduced[i,j])[j]+Y[i]+Y_reduced[i] }     *(M[i,k]+M_reduced[i,k]) * u[k]_t=-grad((A[i,j,k,l]+A_reduced[i,j,k,l]) * grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k]) * u[k])[j]+(C[i,k,l]+C_reduced[i,k,l]) * grad(u[k])[l]+(D[i,k]+D_reduced[i,k] * u[k]-grad(X[i,j]+X_reduced[i,j])[j]+Y[i]+Y_reduced[i]*
3416    
3417     M{A} and M{A_reduced} are of rank four, M{B}, M{B_reduced}, M{C} and     *A* and *A_reduced* are of rank four, *B*, *B_reduced*, *C* and
3418     M{C_reduced} are each of rank three, M{M}, M{M_reduced}, M{D}, M{D_reduced},     *C_reduced* are each of rank three, *M*, *M_reduced*, *D*, *D_reduced*,
3419     M{X_reduced} and M{X} are each of rank two and M{Y} and M{Y_reduced} are of     *X_reduced* and *X* are each of rank two and *Y* and *Y_reduced* are of
3420     rank one. The natural boundary conditions take the form:     rank one. The natural boundary conditions take the form:
3421    
3422     M{n[j]*((A[i,j,k,l]+A_reduced[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k]+X[i,j]+X_reduced[i,j])+(d[i,k]+d_reduced[i,k])*u[k]+y[i]+y_reduced[i]= (m[i,k]+m_reduced[i,k])*u[k]_t}     *n[j]*((A[i,j,k,l]+A_reduced[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k]+X[i,j]+X_reduced[i,j])+(d[i,k]+d_reduced[i,k])*u[k]+y[i]+y_reduced[i]= (m[i,k]+m_reduced[i,k])*u[k]_t*
3423    
3424     The coefficient M{d} and M{m} are of rank two and M{y} is of rank one with     The coefficient *d* and *m* are of rank two and *y* is of rank one with
3425     all in L{FunctionOnBoundary<escript.FunctionOnBoundary>}. The coefficients     all in `FunctionOnBoundary`. The coefficients
3426     M{d_reduced} and M{m_reduced} are of rank two and M{y_reduced} is of rank     *d_reduced* and *m_reduced* are of rank two and *y_reduced* is of rank
3427     one all in L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.     one all in `ReducedFunctionOnBoundary`.
3428    
3429     Constraints take the form     Constraints take the form
3430    
3431     M{u[i]_t=r[i]} where M{q[i]>0}     *u[i]_t=r[i]* where *q[i]>0*
3432    
3433     M{r} and M{q} are each rank one. Notice that at some locations not     *r* and *q* are each rank one. Notice that at some locations not
3434     necessarily all components must have a constraint.     necessarily all components must have a constraint.
3435    
3436     The transport problem is symmetrical if     The transport problem is symmetrical if
3437    
3438        - M{M[i,k]=M[i,k]}        - *M[i,k]=M[i,k]*
3439        - M{M_reduced[i,k]=M_reduced[i,k]}        - *M_reduced[i,k]=M_reduced[i,k]*
3440        - M{A[i,j,k,l]=A[k,l,i,j]}        - *A[i,j,k,l]=A[k,l,i,j]*
3441        - M{A_reduced[i,j,k,l]=A_reduced[k,l,i,j]}        - *A_reduced[i,j,k,l]=A_reduced[k,l,i,j]*
3442        - M{B[i,j,k]=C[k,i,j]}        - *B[i,j,k]=C[k,i,j]*
3443        - M{B_reduced[i,j,k]=C_reduced[k,i,j]}        - *B_reduced[i,j,k]=C_reduced[k,i,j]*
3444        - M{D[i,k]=D[i,k]}        - *D[i,k]=D[i,k]*
3445        - M{D_reduced[i,k]=D_reduced[i,k]}        - *D_reduced[i,k]=D_reduced[i,k]*
3446        - M{m[i,k]=m[k,i]}        - *m[i,k]=m[k,i]*
3447        - M{m_reduced[i,k]=m_reduced[k,i]}        - *m_reduced[i,k]=m_reduced[k,i]*
3448        - M{d[i,k]=d[k,i]}        - *d[i,k]=d[k,i]*
3449        - M{d_reduced[i,k]=d_reduced[k,i]}        - *d_reduced[i,k]=d_reduced[k,i]*
3450    
3451     L{TransportPDE} also supports solution discontinuities over a contact region     `TransportPDE` also supports solution discontinuities over a contact region
3452     in the domain. To specify the conditions across the discontinuity we are     in the domain. To specify the conditions across the discontinuity we are
3453     using the generalised flux M{J} which, in the case of a system of PDEs and     using the generalised flux *J* which, in the case of a system of PDEs and
3454     several components of the solution, is defined as     several components of the solution, is defined as
3455    
3456     M{J[i,j]=(A[i,j,k,l]+A_reduced[[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k]+X[i,j]+X_reduced[i,j]}     *J[i,j]=(A[i,j,k,l]+A_reduced[[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k]+X[i,j]+X_reduced[i,j]*
3457    
3458     For the case of single solution component and single PDE M{J} is defined as     For the case of single solution component and single PDE *J* is defined as
3459    
3460     M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[j]+(B[i]+B_reduced[i])*u+X[i]+X_reduced[i]}     *J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[j]+(B[i]+B_reduced[i])*u+X[i]+X_reduced[i]*
3461    
3462     In the context of discontinuities M{n} denotes the normal on the     In the context of discontinuities *n* denotes the normal on the
3463     discontinuity pointing from side 0 towards side 1 calculated from     discontinuity pointing from side 0 towards side 1 calculated from
3464     L{getNormal<escript.FunctionSpace.getNormal>} of L{FunctionOnContactZero<escript.FunctionOnContactZero>}.     `FunctionSpace.getNormal` of `FunctionOnContactZero`.
3465     For a system of transport problems the contact condition takes the form     For a system of transport problems the contact condition takes the form
3466    
3467     M{n[j]*J0[i,j]=n[j]*J1[i,j]=(y_contact[i]+y_contact_reduced[i])- (d_contact[i,k]+d_contact_reduced[i,k])*jump(u)[k]}     *n[j]*J0[i,j]=n[j]*J1[i,j]=(y_contact[i]+y_contact_reduced[i])- (d_contact[i,k]+d_contact_reduced[i,k])*jump(u)[k]*
3468    
3469     where M{J0} and M{J1} are the fluxes on side 0 and side 1 of the     where *J0* and *J1* are the fluxes on side 0 and side 1 of the
3470     discontinuity, respectively. M{jump(u)}, which is the difference of the     discontinuity, respectively. *jump(u)*, which is the difference of the
3471     solution at side 1 and at side 0, denotes the jump of M{u} across     solution at side 1 and at side 0, denotes the jump of *u* across
3472     discontinuity along the normal calculated by L{jump<util.jump>}.     discontinuity along the normal calculated by `jump`.
3473     The coefficient M{d_contact} is of rank two and M{y_contact} is of rank one     The coefficient *d_contact* is of rank two and *y_contact* is of rank one
3474     both in L{FunctionOnContactZero<escript.FunctionOnContactZero>} or L{FunctionOnContactOne<escript.FunctionOnContactOne>}.     both in `FunctionOnContactZero` or `FunctionOnContactOne`.
3475     The coefficient M{d_contact_reduced} is of rank two and M{y_contact_reduced}     The coefficient *d_contact_reduced* is of rank two and *y_contact_reduced*
3476     is of rank one both in L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}.     is of rank one both in `ReducedFunctionOnContactZero` or `ReducedFunctionOnContactOne`.
3477     In case of a single PDE and a single component solution the contact     In case of a single PDE and a single component solution the contact
3478     condition takes the form     condition takes the form
3479    
3480     M{n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)}     *n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)*
3481    
3482     In this case the coefficient M{d_contact} and M{y_contact} are each scalar     In this case the coefficient *d_contact* and *y_contact* are each scalar
3483     both in L{FunctionOnContactZero<escript.FunctionOnContactZero>} or     both in `FunctionOnContactZero` or
3484     L{FunctionOnContactOne<escript.FunctionOnContactOne>} and the coefficient     `FunctionOnContactOne` and the coefficient
3485     M{d_contact_reduced} and M{y_contact_reduced} are each scalar both in     *d_contact_reduced* and *y_contact_reduced* are each scalar both in
3486     L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or     `ReducedFunctionOnContactZero` or
3487     L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}.     `ReducedFunctionOnContactOne`.
3488    
3489     Typical usage::     Typical usage::
3490    
# Line 2533  class TransportPDE(LinearProblem): Line 3501  class TransportPDE(LinearProblem):
3501       """       """
3502       Initializes a transport problem.       Initializes a transport problem.
3503    
3504       @param domain: domain of the PDE       :param domain: domain of the PDE
3505       @type domain: L{Domain<escript.Domain>}       :type domain: `Domain`
3506       @param numEquations: number of equations. If C{None} the number of       :param numEquations: number of equations. If ``None`` the number of
3507                            equations is extracted from the coefficients.                            equations is extracted from the coefficients.
3508       @param numSolutions: number of solution components. If C{None} the number       :param numSolutions: number of solution components. If ``None`` the number
3509                            of solution components is extracted from the                            of solution components is extracted from the
3510                            coefficients.                            coefficients.
3511       @param debug: if True debug information is printed       :param debug: if True debug information is printed
3512       @param useBackwardEuler: if set the backward Euler scheme is used. Otherwise the Crank-Nicholson scheme is applied. Note that backward Euler scheme will return a safe time step size which is practically infinity as the scheme is unconditional unstable. The Crank-Nicholson scheme provides a higher accuracy but requires to limit the time step size to be stable.       :param useBackwardEuler: if set the backward Euler scheme is used. Otherwise the Crank-Nicolson scheme is applied. Note that backward Euler scheme will return a safe time step size which is practically infinity as the scheme is unconditional unstable. The Crank-Nicolson scheme provides a higher accuracy but requires to limit the time step size to be stable.
3513       @type useBackwardEuler: C{bool}       :type useBackwardEuler: ``bool``
3514       """       """
3515       if useBackwardEuler:       if useBackwardEuler:
3516           self.__useBackwardEuler=True           self.__useBackwardEuler=True
# Line 2586  class TransportPDE(LinearProblem): Line 3554  class TransportPDE(LinearProblem):
3554       """       """
3555       Returns the string representation of the transport problem.       Returns the string representation of the transport problem.
3556    
3557       @return: a simple representation of the transport problem       :return: a simple representation of the transport problem
3558       @rtype: C{str}       :rtype: ``str``
3559       """       """
3560       return "<TransportPDE %d>"%id(self)       return "<TransportPDE %d>"%id(self)
3561    
3562     def useBackwardEuler(self):     def useBackwardEuler(self):
3563        """        """
3564        Returns true if backward Euler is used. Otherwise false is returned.        Returns true if backward Euler is used. Otherwise false is returned.
3565        @rtype: bool        :rtype: bool
3566        """        """
3567        return self.__useBackwardEuler        return self.__useBackwardEuler
3568    
# Line 2603  class TransportPDE(LinearProblem): Line 3571  class TransportPDE(LinearProblem):
3571        """        """
3572        Tests the transport problem for symmetry.        Tests the transport problem for symmetry.
3573    
3574        @param verbose: if set to True or not present a report on coefficients        :param verbose: if set to True or not present a report on coefficients
3575                        which break the symmetry is printed.                        which break the symmetry is printed.
3576        @type verbose: C{bool}        :type verbose: ``bool``
3577        @return:  True if the PDE is symmetric        :return:  True if the PDE is symmetric
3578        @rtype: C{bool}        :rtype: ``bool``
3579        @note: This is a very expensive operation. It should be used for        :note: This is a very expensive operation. It should be used for
3580               degugging only! The symmetry flag is not altered.               degugging only! The symmetry flag is not altered.
3581        """        """
3582        out=True        out=True
# Line 2632  class TransportPDE(LinearProblem): Line 3600  class TransportPDE(LinearProblem):
3600        """        """
3601        Sets new values to coefficients.        Sets new values to coefficients.
3602    
3603        @param coefficients: new values assigned to coefficients        :param coefficients: new values assigned to coefficients
3604        @keyword M: value for coefficient C{M}        :keyword M: value for coefficient ``M``
3605        @type M: any type that can be cast to a L{Data<escript.Data>} object on        :type M: any type that can be cast to a `Data` object on
3606                 L{Function<escript.Function>}                 `Function`
3607        @keyword M_reduced: value for coefficient C{M_reduced}        :keyword M_reduced: value for coefficient ``M_reduced``