/[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 2548 by jfenwick, Mon Jul 20 06:20:06 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-2009 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"""
# Line 21  __url__="https://launchpad.net/escript-f Line 22  __url__="https://launchpad.net/escript-f
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
# Line 52  class SolverOptions(object): Line 53  class SolverOptions(object):
53            
54      Typical usage is      Typical usage is
55            
56      opts=SolverOptions()      ::
57      print opts      
58      opts.resetDiagnostics()        opts=SolverOptions()
59      u=solver(opts)        print opts
60      print "number of iteration steps: =",opts.getDiagnostics("num_iter")        opts.resetDiagnostics()
61              u=solver(opts)
62          print "number of iteration steps: =",opts.getDiagnostics("num_iter")
63      @cvar DEFAULT: The default method used to solve the system of linear equations  
64      @cvar DIRECT: The direct solver based on LDU factorization      :cvar DEFAULT: The default method used to solve the system of linear equations
65      @cvar CHOLEVSKY: The direct solver based on LDLt factorization (can only be applied for symmetric PDEs)      :cvar DIRECT: The direct solver based on LDU factorization
66      @cvar PCG: The preconditioned conjugate gradient method (can only be applied for symmetric PDEs)      :cvar CHOLEVSKY: The direct solver based on LDLt factorization (can only be applied for symmetric PDEs)
67      @cvar CR: The conjugate residual method      :cvar PCG: The preconditioned conjugate gradient method (can only be applied for symmetric PDEs)
68      @cvar CGS: The conjugate gradient square method      :cvar CR: The conjugate residual method
69      @cvar BICGSTAB: The stabilized Bi-Conjugate Gradient method      :cvar CGS: The conjugate gradient square method
70      @cvar TFQMR: Transport Free Quasi Minimal Residual method      :cvar BICGSTAB: The stabilized Bi-Conjugate Gradient method
71      @cvar MINRES: Minimum residual method      :cvar TFQMR: Transpose Free Quasi Minimal Residual method
72      @cvar SSOR: The symmetric over-relaxation method      :cvar MINRES: Minimum residual method
73      @cvar ILU0: The incomplete LU factorization preconditioner with no fill-in      :cvar ILU0: The incomplete LU factorization preconditioner with no fill-in
74      @cvar ILUT: The incomplete LU factorization preconditioner with fill-in      :cvar ILUT: The incomplete LU factorization preconditioner with fill-in
75      @cvar JACOBI: The Jacobi preconditioner      :cvar JACOBI: The Jacobi preconditioner
76      @cvar GMRES: The Gram-Schmidt minimum residual method      :cvar GMRES: The Gram-Schmidt minimum residual method
77      @cvar PRES20: Special GMRES with restart after 20 steps and truncation after 5 residuals      :cvar PRES20: Special GMRES with restart after 20 steps and truncation after 5 residuals
78      @cvar LUMPING: Matrix lumping      :cvar ROWSUM_LUMPING: Matrix lumping using row sum
79      @cvar NO_REORDERING: No matrix reordering allowed      :cvar HRZ_LUMPING: Matrix lumping using the HRZ approach
80      @cvar MINIMUM_FILL_IN: Reorder matrix to reduce fill-in during factorization      :cvar NO_REORDERING: No matrix reordering allowed
81      @cvar NESTED_DISSECTION: Reorder matrix to improve load balancing during factorization      :cvar MINIMUM_FILL_IN: Reorder matrix to reduce fill-in during factorization
82      @cvar PASO: PASO solver package      :cvar NESTED_DISSECTION: Reorder matrix to improve load balancing during factorization
83      @cvar SCSL: SGI SCSL solver library      :cvar PASO: PASO solver package
84      @cvar MKL: Intel's MKL solver library      :cvar SCSL: SGI SCSL solver library
85      @cvar UMFPACK: The UMFPACK library      :cvar MKL: Intel's MKL solver library
86      @cvar TRILINOS: The TRILINOS parallel solver class library from Sandia National Labs      :cvar UMFPACK: The UMFPACK library
87      @cvar ITERATIVE: The default iterative solver      :cvar TRILINOS: The TRILINOS parallel solver class library from Sandia National Labs
88      @cvar AMG: Algebraic Multi Grid      :cvar ITERATIVE: The default iterative solver
89      @cvar REC_ILU: recursive ILU0      :cvar AMG: Algebraic Multi Grid
90      @cvar RILU: relaxed ILU0      :cvar AMLI: Algebraic Multi Level Iteration
91      @cvar GAUSS_SEIDEL: Gauss-Seidel solver      :cvar REC_ILU: recursive ILU0
92      @cvar DEFAULT_REORDERING: the reordering method recommended by the solver      :cvar RILU: relaxed ILU0
93      @cvar SUPER_LU: the Super_LU solver package      :cvar GAUSS_SEIDEL: Gauss-Seidel preconditioner
94      @cvar PASTIX: the Pastix direct solver_package      :cvar DEFAULT_REORDERING: the reordering method recommended by the solver
95      @cvar YAIR_SHAPIRA_COARSENING: AMG coarsening method by Yair-Shapira      :cvar SUPER_LU: the Super_LU solver package
96      @cvar RUGE_STUEBEN_COARSENING: AMG coarsening method by Ruge and Stueben      :cvar PASTIX: the Pastix direct solver_package
97      @cvar AGGREGATION_COARSENING: AMG coarsening using (symmetric) aggregation      :cvar YAIR_SHAPIRA_COARSENING: AMG and AMLI coarsening method by Yair-Shapira
98      @cvar MIN_COARSE_MATRIX_SIZE: minimum size of the coarsest level matrix to use direct solver.      :cvar RUGE_STUEBEN_COARSENING: AMG and AMLI coarsening method by Ruge and Stueben
99      @cvar NO_PRECONDITIONER: no preconditioner is applied.      :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      DEFAULT= 0
105      DIRECT= 1      DIRECT= 1
# Line 104  class SolverOptions(object): Line 108  class SolverOptions(object):
108      CR= 4      CR= 4
109      CGS= 5      CGS= 5
110      BICGSTAB= 6      BICGSTAB= 6
     SSOR= 7  
111      ILU0= 8      ILU0= 8
112      ILUT= 9      ILUT= 9
113      JACOBI= 10      JACOBI= 10
114      GMRES= 11      GMRES= 11
115      PRES20= 12      PRES20= 12
116      LUMPING= 13      LUMPING= 13
117        ROWSUM_LUMPING= 13
118        HRZ_LUMPING= 14
119      NO_REORDERING= 17      NO_REORDERING= 17
120      MINIMUM_FILL_IN= 18      MINIMUM_FILL_IN= 18
121      NESTED_DISSECTION= 19      NESTED_DISSECTION= 19
# Line 134  class SolverOptions(object): Line 139  class SolverOptions(object):
139      AGGREGATION_COARSENING=35      AGGREGATION_COARSENING=35
140      NO_PRECONDITIONER=36      NO_PRECONDITIONER=36
141      MIN_COARSE_MATRIX_SIZE=37      MIN_COARSE_MATRIX_SIZE=37
142            AMLI=38
143        STANDARD_COARSENING=39
144    
145      def __init__(self):      def __init__(self):
146          self.setLevelMax()          self.setLevelMax()
147          self.setCoarseningThreshold()          self.setCoarseningThreshold()
148            self.setSmoother()
149          self.setNumSweeps()          self.setNumSweeps()
150          self.setNumPreSweeps()          self.setNumPreSweeps()
151          self.setNumPostSweeps()          self.setNumPostSweeps()
# Line 161  class SolverOptions(object): Line 169  class SolverOptions(object):
169          self.setCoarsening()          self.setCoarsening()
170          self.setMinCoarseMatrixSize()          self.setMinCoarseMatrixSize()
171          self.setRelaxationFactor()                  self.setRelaxationFactor()        
172            self.setLocalPreconditionerOff()
173          self.resetDiagnostics(all=True)          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):      def __str__(self):
183          return self.getSummary()          return self.getSummary()
# Line 176  class SolverOptions(object): Line 192  class SolverOptions(object):
192          out+="\nAbsolute tolerance = %e"%self.getAbsoluteTolerance()          out+="\nAbsolute tolerance = %e"%self.getAbsoluteTolerance()
193          out+="\nSymmetric problem = %s"%self.isSymmetric()          out+="\nSymmetric problem = %s"%self.isSymmetric()
194          out+="\nMaximum number of iteration steps = %s"%self.getIterMax()          out+="\nMaximum number of iteration steps = %s"%self.getIterMax()
195          out+="\nInner tolerance = %e"%self.getInnerTolerance()          # out+="\nInner tolerance = %e"%self.getInnerTolerance()
196          out+="\nAdapt innner tolerance = %s"%self.adaptInnerTolerance()          # out+="\nAdapt innner tolerance = %s"%self.adaptInnerTolerance()
197            
198          if self.getPackage() == self.PASO:          if self.getPackage() == self.PASO:
199              out+="\nSolver method = %s"%self.getName(self.getSolverMethod())              out+="\nSolver method = %s"%self.getName(self.getSolverMethod())
# Line 186  class SolverOptions(object): Line 202  class SolverOptions(object):
202                  out+="\nRestart  = %s"%self.getRestart()                  out+="\nRestart  = %s"%self.getRestart()
203              if self.getSolverMethod() == self.AMG:              if self.getSolverMethod() == self.AMG:
204                  out+="\nNumber of pre / post sweeps = %s / %s, %s"%(self.getNumPreSweeps(), self.getNumPostSweeps(), self.getNumSweeps())                  out+="\nNumber of pre / post sweeps = %s / %s, %s"%(self.getNumPreSweeps(), self.getNumPostSweeps(), self.getNumSweeps())
205                  out+="\nMaximum number of levels = %s"%self.LevelMax()                  out+="\nMaximum number of levels = %s"%self.getLevelMax()
206                  out+="\nCoarsening threshold = %e"%self.getCoarseningThreshold()                  out+="\nCoarsening threshold = %e"%self.getCoarseningThreshold()
207                  out+="\Coarsening method = %s"%self.getName(self.getCoarsening())                  out+="\nCoarsening method = %s"%self.getName(self.getCoarsening())
208              out+="\nPreconditioner = %s"%self.getName(self.getPreconditioner())              out+="\nPreconditioner = %s"%self.getName(self.getPreconditioner())
209                out+="\nApply preconditioner locally = %s"%self.useLocalPreconditioner()
210              if self.getPreconditioner() == self.AMG:              if self.getPreconditioner() == self.AMG:
211                  out+="\nMaximum number of levels = %s"%self.LevelMax()                  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())                  out+="\nCoarsening method = %s"%self.getName(self.getCoarsening())
225                  out+="\nCoarsening threshold = %e"%self.getMinCoarseMatrixSize()                  out+="\nCoarsening threshold = %e"%self.getMinCoarseMatrixSize()
226                  out+="\nMinimum size of the coarsest level matrix = %e"%self.getCoarseningThreshold()                  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())                  out+="\nNumber of pre / post sweeps = %s / %s, %s"%(self.getNumPreSweeps(), self.getNumPostSweeps(), self.getNumSweeps())
228              if self.getPreconditioner() == self.GAUSS_SEIDEL:          if self.getPreconditioner() == self.GAUSS_SEIDEL:
229                  out+="\nNumber of sweeps = %s"%self.getNumSweeps()                  out+="\nNumber of sweeps = %s"%self.getNumSweeps()
230              if self.getPreconditioner() == self.ILUT:              if self.getPreconditioner() == self.ILUT:
231                  out+="\nDrop tolerance = %e"%self.getDropTolerance()                  out+="\nDrop tolerance = %e"%self.getDropTolerance()
# Line 209  class SolverOptions(object): Line 238  class SolverOptions(object):
238          """          """
239          returns the name of a given key          returns the name of a given key
240                    
241          @param key: a valid key          :param key: a valid key
242          """          """
243          if key == self.DEFAULT: return "DEFAULT"          if key == self.DEFAULT: return "DEFAULT"
244          if key == self.DIRECT: return "DIRECT"          if key == self.DIRECT: return "DIRECT"
# Line 218  class SolverOptions(object): Line 247  class SolverOptions(object):
247          if key == self.CR: return "CR"          if key == self.CR: return "CR"
248          if key == self.CGS: return "CGS"          if key == self.CGS: return "CGS"
249          if key == self.BICGSTAB: return "BICGSTAB"          if key == self.BICGSTAB: return "BICGSTAB"
         if key == self.SSOR: return "SSOR"  
250          if key == self.ILU0: return "ILU0:"          if key == self.ILU0: return "ILU0:"
251          if key == self.ILUT: return "ILUT"          if key == self.ILUT: return "ILUT"
252          if key == self.JACOBI: return "JACOBI"          if key == self.JACOBI: return "JACOBI"
253          if key == self.GMRES: return "GMRES"          if key == self.GMRES: return "GMRES"
254          if key == self.PRES20: return "PRES20"          if key == self.PRES20: return "PRES20"
255          if key == self.LUMPING: return "LUMPING"          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"          if key == self.NO_REORDERING: return "NO_REORDERING"
258          if key == self.MINIMUM_FILL_IN: return "MINIMUM_FILL_IN"          if key == self.MINIMUM_FILL_IN: return "MINIMUM_FILL_IN"
259          if key == self.NESTED_DISSECTION: return "NESTED_DISSECTION"          if key == self.NESTED_DISSECTION: return "NESTED_DISSECTION"
# Line 233  class SolverOptions(object): Line 262  class SolverOptions(object):
262          if key == self.ITERATIVE: return "ITERATIVE"          if key == self.ITERATIVE: return "ITERATIVE"
263          if key == self.PASO: return "PASO"          if key == self.PASO: return "PASO"
264          if key == self.AMG: return "AMG"          if key == self.AMG: return "AMG"
265            if key == self.AMLI: return "AMLI"
266          if key == self.REC_ILU: return "REC_ILU"          if key == self.REC_ILU: return "REC_ILU"
267          if key == self.TRILINOS: return "TRILINOS"          if key == self.TRILINOS: return "TRILINOS"
268          if key == self.NONLINEAR_GMRES: return "NONLINEAR_GMRES"          if key == self.NONLINEAR_GMRES: return "NONLINEAR_GMRES"
# Line 245  class SolverOptions(object): Line 275  class SolverOptions(object):
275          if key == self.PASTIX: return "PASTIX"          if key == self.PASTIX: return "PASTIX"
276          if key == self.YAIR_SHAPIRA_COARSENING: return "YAIR_SHAPIRA_COARSENING"          if key == self.YAIR_SHAPIRA_COARSENING: return "YAIR_SHAPIRA_COARSENING"
277          if key == self.RUGE_STUEBEN_COARSENING: return "RUGE_STUEBEN_COARSENING"          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"          if key == self.AGGREGATION_COARSENING: return "AGGREGATION_COARSENING"
280          if key == self.NO_PRECONDITIONER: return "NO_PRECONDITIONER"          if key == self.NO_PRECONDITIONER: return "NO_PRECONDITIONER"
         if key == self.MIN_COARSE_MATRIX_SIZE: return "MIN_COARSE_MATRIX_SIZE"  
281                    
282      def resetDiagnostics(self,all=False):      def resetDiagnostics(self,all=False):
283          """          """
284          resets the diagnostics          resets the diagnostics
285                    
286          @param all: if C{all} is C{True} all diagnostics including accumulative counters are reset.          :param all: if ``all`` is ``True`` all diagnostics including accumulative counters are reset.
287          @type all: C{bool}          :type all: ``bool``
288          """          """
289          self.__num_iter=None          self.__num_iter=None
290          self.__num_level=None          self.__num_level=None
291          self.__num_inner_iter=None          self.__num_inner_iter=None
292          self.__time=None          self.__time=None
293          self.__set_up_time=None          self.__set_up_time=None
294            self.__net_time=None
295          self.__residual_norm=None          self.__residual_norm=None
296          self.__converged=None          self.__converged=None
297            self.__preconditioner_size=-1
298            self.__time_step_backtracking_used=None
299          if all:          if all:
300              self.__cum_num_inner_iter=0              self.__cum_num_inner_iter=0
301              self.__cum_num_iter=0              self.__cum_num_iter=0
302              self.__cum_time=0              self.__cum_time=0
303              self.__cum_set_up_time=0              self.__cum_set_up_time=0
304                self.__cum_net_time=0
305    
306      def _updateDiagnostics(self, name, value):      def _updateDiagnostics(self, name, value):
307          """          """
308          Updates diagnostic information          Updates diagnostic information
309                    
310          @param name: name of  diagnostic information          :param name: name of  diagnostic information
311          @type name: C{str} in the list "num_iter", "num_level", "num_inner_iter", "time", "set_up_time", "residual_norm", "converged".          :type name: ``str`` in the list "num_iter", "num_level", "num_inner_iter", "time", "set_up_time", "net_time", "residual_norm", "converged".
312          @param vale: new value of the diagnostic information          :param value: new value of the diagnostic information
313          @note: this function is used by a solver to report diagnostics informations.          :note: this function is used by a solver to report diagnostics informations.
314          """          """
315          if name == "num_iter":          if name == "num_iter":
316              self.__num_iter=int(value)              self.__num_iter=int(value)
317              self.__cum_num_iter+=self.__num_iter              self.__cum_num_iter+=self.__num_iter
318          if name == "num_level":          if name == "num_level":
319              self.__num_iter=int(value)              self.__num_level=int(value)
320          if name == "num_inner_iter":          if name == "num_inner_iter":
321              self.__num_inner_iter=int(value)              self.__num_inner_iter=int(value)
322              self.__cum_num_inner_iter+=self.__num_inner_iter              self.__cum_num_inner_iter+=self.__num_inner_iter
# Line 292  class SolverOptions(object): Line 326  class SolverOptions(object):
326          if name == "set_up_time":          if name == "set_up_time":
327              self.__set_up_time=float(value)              self.__set_up_time=float(value)
328              self.__cum_set_up_time+=self.__set_up_time              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":          if name == "residual_norm":
333              self.__residual_norm=float(value)              self.__residual_norm=float(value)
334          if name == "converged":          if name == "converged":
335              self.__converged = (value == True)              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):      def getDiagnostics(self, name):
343          """          """
344          Returns the diagnostic information C{name}          Returns the diagnostic information ``name``. Possible values are:
345                    
346          @param name: name of diagnostic information where      - "num_iter": the number of iteration steps
         - "num_iter": the number of iteration steps  
347          - "cum_num_iter": the cumulative number of iteration steps          - "cum_num_iter": the cumulative number of iteration steps
348          - "num_level": the number of level in multi level solver          - "num_level": the number of level in multi level solver
349          - "num_inner_iter": the number of inner iteration steps          - "num_inner_iter": the number of inner iteration steps
# Line 310  class SolverOptions(object): Line 352  class SolverOptions(object):
352          - "cum_time": cumulative execution time          - "cum_time": cumulative execution time
353          - "set_up_time": time to set up of the solver, typically this includes factorization and reordering          - "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          - "cum_set_up_time": cumulative time to set up of the solver
355          - "residual_norm": norm of the final residual          - "net_time": net execution time, excluding setup time for the solver and execution time for preconditioner
356          - "converged": return self.__converged              - "cum_net_time": cumulative net execution time
357          @type name: C{str} in the list "num_iter", "num_level", "num_inner_iter", "time", "set_up_time", "residual_norm", "converged".          - "preconditioner_size": size of preconditioner [Bytes]
358          @return: requested value. C{None} is returned if the value is undefined.          - "converged": return True if solution has converged.
359          @note: If the solver has thrown an exception diagnostic values have an undefined status.          - "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          if name == "num_iter": return self.__num_iter
370          elif name == "cum_num_iter": return self.__cum_num_iter          elif name == "cum_num_iter": return self.__cum_num_iter
# Line 325  class SolverOptions(object): Line 375  class SolverOptions(object):
375          elif name == "cum_time": return self.__cum_time          elif name == "cum_time": return self.__cum_time
376          elif name == "set_up_time": return self.__set_up_time          elif name == "set_up_time": return self.__set_up_time
377          elif name == "cum_set_up_time": return self.__cum_set_up_time          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          elif name == "residual_norm": return self.__residual_norm
381          elif name == "converged": return self.__converged                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:          else:
387              raise ValueError,"unknown diagnostic item %s"%name              raise ValueError,"unknown diagnostic item %s"%name
388      def hasConverged(self):      def hasConverged(self):
389          """          """
390          Returns C{True} if the last solver call has been finalized successfully.          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.          :note: if an exception has been thrown by the solver the status of this flag is undefined.
392          """          """
393          return self.getDiagnostics("converged")          return self.getDiagnostics("converged")
394      def setCoarsening(self,method=0):      def setCoarsening(self,method=0):
395          """          """
396          Sets the key of the coarsening method to be applied in AMG.          Sets the key of the coarsening method to be applied in AMG or AMLI
397    
398          @param method: selects the coarsening method .          :param method: selects the coarsening method .
399          @type method: in {SolverOptions.DEFAULT}, L{SolverOptions.YAIR_SHAPIRA_COARSENING},          :type method: in {SolverOptions.DEFAULT}, `SolverOptions.YAIR_SHAPIRA_COARSENING`,  `SolverOptions.RUGE_STUEBEN_COARSENING`, `SolverOptions.AGGREGATION_COARSENING`
         L{SolverOptions.RUGE_STUEBEN_COARSENING}, L{SolverOptions.AGGREGATION_COARSENING}  
400          """          """
401      if method==None: method=0      if method==None: method=0
402          if not method in [self.DEFAULT, self.YAIR_SHAPIRA_COARSENING, self.RUGE_STUEBEN_COARSENING, self.AGGREGATION_COARSENING]:          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               raise ValueError,"unknown coarsening method %s"%method
404          self.__coarsening=method          self.__coarsening=method
405            
406      def getCoarsening(self):      def getCoarsening(self):
407          """          """
408          Returns the key of the coarsening algorithm to be applied AMG.          Returns the key of the coarsening algorithm to be applied AMG or AMLI
409    
410          @rtype: in the list L{SolverOptions.DEFAULT}, L{SolverOptions.YAIR_SHAPIRA_COARSENING},          :rtype: in the list `SolverOptions.DEFAULT`, `SolverOptions.YAIR_SHAPIRA_COARSENING`, `SolverOptions.RUGE_STUEBEN_COARSENING`, `SolverOptions.AGGREGATION_COARSENING`
         L{SolverOptions.RUGE_STUEBEN_COARSENING}, L{SolverOptions.AGGREGATION_COARSENING}  
411          """          """
412          return self.__coarsening          return self.__coarsening
413                
414      def setMinCoarseMatrixSize(self,size=500):      def setMinCoarseMatrixSize(self,size=500):
415          """          """
416          Sets the minumum size of the coarsest level matrix in AMG.          Sets the minumum size of the coarsest level matrix in AMG or AMLI
417    
418          @param size: minumum size of the coarsest level matrix .          :param size: minumum size of the coarsest level matrix .
419          @type size: positive C{int} or C{None}          :type size: positive ``int`` or ``None``
420          """          """
421        if size==None: size=500
422          size=int(size)          size=int(size)
423          if size<0:          if size<0:
424             raise ValueError,"minumum size of the coarsest level matrix must be non-negative."             raise ValueError,"minumum size of the coarsest level matrix must be non-negative."
     if size==None: size=500  
425          self.__MinCoarseMatrixSize=size          self.__MinCoarseMatrixSize=size
426                    
427      def getMinCoarseMatrixSize(self):      def getMinCoarseMatrixSize(self):
428          """          """
429          Returns the minumum size of the coarsest level matrix in AMG.          Returns the minumum size of the coarsest level matrix in AMG or AMLI
430    
431          @rtype: C{int}          :rtype: ``int``
432          """          """
433          return self.__MinCoarseMatrixSize          return self.__MinCoarseMatrixSize
434                
# Line 382  class SolverOptions(object): Line 436  class SolverOptions(object):
436          """          """
437          Sets the preconditioner to be used.          Sets the preconditioner to be used.
438    
439          @param preconditioner: key of the preconditioner to be used.          :param preconditioner: key of the preconditioner to be used.
440          @type preconditioner: in L{SolverOptions.SSOR}, L{SolverOptions.ILU0}, L{SolverOptions.ILUT}, L{SolverOptions.JACOBI},          :type preconditioner: in `SolverOptions.ILU0`, `SolverOptions.ILUT`, `SolverOptions.JACOBI`,
441                                      L{SolverOptions.AMG}, L{SolverOptions.REC_ILU}, L{SolverOptions.GAUSS_SEIDEL}, L{SolverOptions.RILU},                                      `SolverOptions.AMG`, `SolverOptions.AMLI`, `SolverOptions.REC_ILU`, `SolverOptions.GAUSS_SEIDEL`, `SolverOptions.RILU`,
442                                      L{SolverOptions.NO_PRECONDITIONER}                                      `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          :note: Not all packages support all preconditioner. It can be assumed that a package makes a reasonable choice if it encounters an unknown preconditioner.
         an unknown preconditioner.  
444          """          """
445      if preconditioner==None: preconditioner=10      if preconditioner==None: preconditioner=10
446          if not preconditioner in [ SolverOptions.SSOR, SolverOptions.ILU0, SolverOptions.ILUT, SolverOptions.JACOBI,          if not preconditioner in [  SolverOptions.ILU0, SolverOptions.ILUT, SolverOptions.JACOBI,
447                                      SolverOptions.AMG, SolverOptions.REC_ILU, SolverOptions.GAUSS_SEIDEL, SolverOptions.RILU,                                      SolverOptions.AMG, SolverOptions.AMLI, SolverOptions.REC_ILU, SolverOptions.GAUSS_SEIDEL, SolverOptions.RILU,
448                                      SolverOptions.NO_PRECONDITIONER] :                                      SolverOptions.NO_PRECONDITIONER] :
449               raise ValueError,"unknown preconditioner %s"%preconditioner               raise ValueError,"unknown preconditioner %s"%preconditioner
450          self.__preconditioner=preconditioner              self.__preconditioner=preconditioner    
# Line 399  class SolverOptions(object): Line 452  class SolverOptions(object):
452          """          """
453          Returns key of the preconditioner to be used.          Returns key of the preconditioner to be used.
454    
455          @rtype: in the list L{SolverOptions.SSOR}, L{SolverOptions.ILU0}, L{SolverOptions.ILUT}, L{SolverOptions.JACOBI},          :rtype: in the list `SolverOptions.ILU0`, `SolverOptions.ILUT`, `SolverOptions.JACOBI`, SolverOptions.AMLI,
456                                      L{SolverOptions.AMG}, L{SolverOptions.REC_ILU}, L{SolverOptions.GAUSS_SEIDEL}, L{SolverOptions.RILU},                                      `SolverOptions.AMG`, `SolverOptions.REC_ILU`, `SolverOptions.GAUSS_SEIDEL`, `SolverOptions.RILU`,
457                                      L{SolverOptions.NO_PRECONDITIONER}                                      `SolverOptions.NO_PRECONDITIONER`
458          """          """
459          return self.__preconditioner          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):      def setSolverMethod(self, method=0):
480          """          """
481          Sets the solver method to be used. Use C{method}=C{SolverOptions.DIRECT} to indicate that a direct rather than an iterative          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 C{method}=C{SolverOptions.ITERATIVE} to indicate that an iterative rather than a direct          solver should be used and Use ``method``=``SolverOptions.ITERATIVE`` to indicate that an iterative rather than a direct
483          solver should be used.          solver should be used.
484    
485          @param method: key of the solver method to be used.          :param method: key of the solver method to be used.
486          @type method: in L{SolverOptions.DEFAULT}, L{SolverOptions.DIRECT}, L{SolverOptions.CHOLEVSKY}, L{SolverOptions.PCG},          :type method: in `SolverOptions.DEFAULT`, `SolverOptions.DIRECT`, `SolverOptions.CHOLEVSKY`, `SolverOptions.PCG`,
487                          L{SolverOptions.CR}, L{SolverOptions.CGS}, L{SolverOptions.BICGSTAB}, L{SolverOptions.SSOR},                          `SolverOptions.CR`, `SolverOptions.CGS`, `SolverOptions.BICGSTAB`,
488                          L{SolverOptions.GMRES}, L{SolverOptions.PRES20}, L{SolverOptions.LUMPING}, L{SolverOptions.ITERATIVE},                          `SolverOptions.GMRES`, `SolverOptions.PRES20`, `SolverOptions.ROWSUM_LUMPING`, `SolverOptions.HRZ_LUMPING`, `SolverOptions.ITERATIVE`,
489                          L{SolverOptions.AMG}, L{SolverOptions.NONLINEAR_GMRES}, L{SolverOptions.TFQMR}, L{SolverOptions.MINRES},                          `SolverOptions.NONLINEAR_GMRES`, `SolverOptions.TFQMR`, `SolverOptions.MINRES`
490                          L{SolverOptions.GAUSS_SEIDEL}          :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.
         @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      if method==None: method=0
493          if not method in [ SolverOptions.DEFAULT, SolverOptions.DIRECT, SolverOptions.CHOLEVSKY, SolverOptions.PCG,          if not method in [ SolverOptions.DEFAULT, SolverOptions.DIRECT, SolverOptions.CHOLEVSKY, SolverOptions.PCG,
494                             SolverOptions.CR, SolverOptions.CGS, SolverOptions.BICGSTAB, SolverOptions.SSOR,                             SolverOptions.CR, SolverOptions.CGS, SolverOptions.BICGSTAB,
495                             SolverOptions.GMRES, SolverOptions.PRES20, SolverOptions.LUMPING, SolverOptions.ITERATIVE, SolverOptions.AMG,                             SolverOptions.GMRES, SolverOptions.PRES20, SolverOptions.ROWSUM_LUMPING, SolverOptions.HRZ_LUMPING,
496                             SolverOptions.NONLINEAR_GMRES, SolverOptions.TFQMR, SolverOptions.MINRES, SolverOptions.GAUSS_SEIDEL]:                             SolverOptions.ITERATIVE,
497                               SolverOptions.NONLINEAR_GMRES, SolverOptions.TFQMR, SolverOptions.MINRES ]:
498               raise ValueError,"unknown solver method %s"%method               raise ValueError,"unknown solver method %s"%method
499          self.__method=method          self.__method=method
500      def getSolverMethod(self):      def getSolverMethod(self):
501          """          """
502          Returns key of the solver method to be used.          Returns key of the solver method to be used.
503    
504          @rtype: in the list L{SolverOptions.DEFAULT}, L{SolverOptions.DIRECT}, L{SolverOptions.CHOLEVSKY}, L{SolverOptions.PCG},          :rtype: in the list `SolverOptions.DEFAULT`, `SolverOptions.DIRECT`, `SolverOptions.CHOLEVSKY`, `SolverOptions.PCG`,
505                          L{SolverOptions.CR}, L{SolverOptions.CGS}, L{SolverOptions.BICGSTAB}, L{SolverOptions.SSOR},                          `SolverOptions.CR`, `SolverOptions.CGS`, `SolverOptions.BICGSTAB`,
506                          L{SolverOptions.GMRES}, L{SolverOptions.PRES20}, L{SolverOptions.LUMPING}, L{SolverOptions.ITERATIVE},                          `SolverOptions.GMRES`, `SolverOptions.PRES20`, `SolverOptions.ROWSUM_LUMPING`, `SolverOptions.HRZ_LUMPING`, `SolverOptions.ITERATIVE`,
507                          L{SolverOptions.AMG}, L{SolverOptions.NONLINEAR_GMRES}, L{SolverOptions.TFQMR}, L{SolverOptions.MINRES},                          `SolverOptions.NONLINEAR_GMRES`, `SolverOptions.TFQMR`, `SolverOptions.MINRES`
                         L{SolverOptions.GAUSS_SEIDEL}  
508          """          """
509          return self.__method          return self.__method
510                    
# Line 442  class SolverOptions(object): Line 512  class SolverOptions(object):
512          """          """
513          Sets the solver package to be used as a solver.            Sets the solver package to be used as a solver.  
514    
515          @param package: key of the solver package to be used.          :param package: key of the solver package to be used.
516          @type package: in L{SolverOptions.DEFAULT}, L{SolverOptions.PASO}, L{SolverOptions.SUPER_LU}, L{SolverOptions.PASTIX}, L{SolverOptions.MKL}, L{SolverOptions.UMFPACK}, L{SolverOptions.TRILINOS}          :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.          :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      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]:          if not package in [SolverOptions.DEFAULT, SolverOptions.PASO, SolverOptions.SUPER_LU, SolverOptions.PASTIX, SolverOptions.MKL, SolverOptions.UMFPACK, SolverOptions.TRILINOS]:
# Line 454  class SolverOptions(object): Line 524  class SolverOptions(object):
524          """          """
525          Returns the solver package key          Returns the solver package key
526    
527          @rtype: in the list L{SolverOptions.DEFAULT}, L{SolverOptions.PASO}, L{SolverOptions.SUPER_LU}, L{SolverOptions.PASTIX}, L{SolverOptions.MKL}, L{SolverOptions.UMFPACK}, L{SolverOptions.TRILINOS}          :rtype: in the list `SolverOptions.DEFAULT`, `SolverOptions.PASO`, `SolverOptions.SUPER_LU`, `SolverOptions.PASTIX`, `SolverOptions.MKL`, `SolverOptions.UMFPACK`, `SolverOptions.TRILINOS`
528          """          """
529          return self.__package          return self.__package
530      def setReordering(self,ordering=30):      def setReordering(self,ordering=30):
# Line 462  class SolverOptions(object): Line 532  class SolverOptions(object):
532          Sets the key of the reordering method to be applied if supported by the solver. Some direct solvers support reordering          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.          to optimize compute time and storage use during elimination.
534    
535          @param ordering: selects the reordering strategy.          :param ordering: selects the reordering strategy.
536          @type ordering: in L{SolverOptions.NO_REORDERING}, L{SolverOptions.NO_REORDERING},          :type ordering: in 'SolverOptions.NO_REORDERING', 'SolverOptions.MINIMUM_FILL_IN', 'SolverOptions.NESTED_DISSECTION', 'SolverOptions.DEFAULT_REORDERING'
         L{SolverOptions.NO_REORDERING}, L{SolverOptions.DEFAULT_REORDERING}  
537          """          """
538          if not ordering in [self.NO_REORDERING, self.MINIMUM_FILL_IN, self.NESTED_DISSECTION, self.DEFAULT_REORDERING]:          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               raise ValueError,"unknown reordering strategy %s"%ordering
# Line 473  class SolverOptions(object): Line 542  class SolverOptions(object):
542          """          """
543          Returns the key of the reordering method to be applied if supported by the solver.          Returns the key of the reordering method to be applied if supported by the solver.
544    
545          @rtype: in the list L{SolverOptions.NO_REORDERING}, L{SolverOptions.NO_REORDERING},          :rtype: ordering in 'SolverOptions.NO_REORDERING', 'SolverOptions.MINIMUM_FILL_IN', 'SolverOptions.NESTED_DISSECTION', 'SolverOptions.DEFAULT_REORDERING'
         L{SolverOptions.NO_REORDERING}, L{SolverOptions.DEFAULT_REORDERING}  
546          """          """
547          return self.__reordering          return self.__reordering
548      def setRestart(self,restart=None):      def setRestart(self,restart=None):
549          """          """
550          Sets the number of iterations steps after which GMRES is performing a restart.          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 C{None} no          :param restart: number of iteration steps after which to perform a restart. If equal to ``None`` no restart is performed.
553                          restart is performed.          :type restart: ``int`` or ``None``
         @type restart: C{int} or C{None}  
554          """          """
555          if restart == None:          if restart == None:
556              self.__restart=restart              self.__restart=restart
# Line 496  class SolverOptions(object): Line 563  class SolverOptions(object):
563      def getRestart(self):      def getRestart(self):
564          """          """
565          Returns the number of iterations steps after which GMRES is performing a restart.          Returns the number of iterations steps after which GMRES is performing a restart.
566          If C{None} is returned no restart is performed.          If ``None`` is returned no restart is performed.
567    
568          @rtype: C{int} or C{None}          :rtype: ``int`` or ``None``
569          """          """
570          if self.__restart < 0:          if self.__restart < 0:
571              return None              return None
# Line 510  class SolverOptions(object): Line 577  class SolverOptions(object):
577              return -1              return -1
578              else:              else:
579              return r              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):      def setTruncation(self,truncation=20):
602          """          """
603          Sets the number of residuals in GMRES to be stored for orthogonalization.  The more residuals are stored          Sets the number of residuals in GMRES to be stored for orthogonalization.  The more residuals are stored
604          the faster GMRES converged but          the faster GMRES converged but
605    
606          @param truncation: truncation          :param truncation: truncation
607          @type truncation: C{int}          :type truncation: ``int``
608          """          """
609          truncation=int(truncation)          truncation=int(truncation)
610          if truncation<1:          if truncation<1:
# Line 526  class SolverOptions(object): Line 614  class SolverOptions(object):
614          """          """
615          Returns the number of residuals in GMRES to be stored for orthogonalization          Returns the number of residuals in GMRES to be stored for orthogonalization
616    
617          @rtype: C{int}          :rtype: ``int``
618          """          """
619          return self.__truncation          return self.__truncation
620      def setInnerIterMax(self,iter_max=10):      def setInnerIterMax(self,iter_max=10):
621          """          """
622          Sets the maximum number of iteration steps for the inner iteration.          Sets the maximum number of iteration steps for the inner iteration.
623    
624          @param iter_max: maximum number of inner iterations          :param iter_max: maximum number of inner iterations
625          @type iter_max: C{int}          :type iter_max: ``int``
626          """          """
627          iter_max=int(iter_max)          iter_max=int(iter_max)
628          if iter_max<1:          if iter_max<1:
# Line 544  class SolverOptions(object): Line 632  class SolverOptions(object):
632          """          """
633          Returns maximum number of inner iteration steps          Returns maximum number of inner iteration steps
634    
635          @rtype: C{int}          :rtype: ``int``
636          """          """
637          return self.__inner_iter_max          return self.__inner_iter_max
638      def setIterMax(self,iter_max=100000):      def setIterMax(self,iter_max=100000):
639          """          """
640          Sets the maximum number of iteration steps          Sets the maximum number of iteration steps
641    
642          @param iter_max: maximum number of iteration steps          :param iter_max: maximum number of iteration steps
643          @type iter_max: C{int}          :type iter_max: ``int``
644          """          """
645          iter_max=int(iter_max)          iter_max=int(iter_max)
646          if iter_max<1:          if iter_max<1:
# Line 562  class SolverOptions(object): Line 650  class SolverOptions(object):
650          """          """
651          Returns maximum number of iteration steps          Returns maximum number of iteration steps
652    
653          @rtype: C{int}          :rtype: ``int``
654          """          """
655          return self.__iter_max          return self.__iter_max
656      def setLevelMax(self,level_max=10):      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          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          :param level_max: maximum number of levels
661          @type level_max: C{int}          :type level_max: ``int``
662          """          """
663          level_max=int(level_max)          level_max=int(level_max)
664          if level_max<0:          if level_max<0:
# Line 580  class SolverOptions(object): Line 668  class SolverOptions(object):
668          """          """
669          Returns the maximum number of coarsening levels to be used in an algebraic multi level solver or preconditioner          Returns the maximum number of coarsening levels to be used in an algebraic multi level solver or preconditioner
670    
671          @rtype: C{int}          :rtype: ``int``
672          """          """
673          return self.__level_max          return self.__level_max
674      def setCoarseningThreshold(self,theta=0.05):      def setCoarseningThreshold(self,theta=0.25):
675          """          """
676          Sets the threshold for coarsening in the algebraic multi level solver or preconditioner          Sets the threshold for coarsening in the algebraic multi level solver or preconditioner
677    
678          @param theta: threshold for coarsening          :param theta: threshold for coarsening
679          @type theta: positive C{float}          :type theta: positive ``float``
680          """          """
681          theta=float(theta)          theta=float(theta)
682          if theta<0 or theta>1:          if theta<0 or theta>1:
# Line 598  class SolverOptions(object): Line 686  class SolverOptions(object):
686          """          """
687          Returns the threshold for coarsening in the algebraic multi level solver or preconditioner          Returns the threshold for coarsening in the algebraic multi level solver or preconditioner
688    
689          @rtype: C{float}          :rtype: ``float``
690          """          """
691          return self.__coarsening_threshold          return self.__coarsening_threshold
692      def setNumSweeps(self,sweeps=2):      def setNumSweeps(self,sweeps=1):
693          """          """
694          Sets the number of sweeps in a Jacobi or Gauss-Seidel/SOR preconditioner.          Sets the number of sweeps in a Jacobi or Gauss-Seidel/SOR preconditioner.
695    
696          @param sweeps: number of sweeps          :param sweeps: number of sweeps
697          @type theta: positive C{int}          :type sweeps: positive ``int``
698          """          """
699          sweeps=int(sweeps)          sweeps=int(sweeps)
700          if sweeps<1:          if sweeps<1:
# Line 616  class SolverOptions(object): Line 704  class SolverOptions(object):
704          """          """
705          Returns the number of sweeps in a Jacobi or Gauss-Seidel/SOR preconditioner.          Returns the number of sweeps in a Jacobi or Gauss-Seidel/SOR preconditioner.
706    
707          @rtype: C{int}          :rtype: ``int``
708          """          """
709          return self.__sweeps          return self.__sweeps
710      def setNumPreSweeps(self,sweeps=2):      def setNumPreSweeps(self,sweeps=2):
711          """          """
712          Sets the number of sweeps in the pre-smoothing step of a multi level solver or preconditioner          Sets the number of sweeps in the pre-smoothing step of a multi level solver or preconditioner
713    
714          @param sweeps: number of sweeps          :param sweeps: number of sweeps
715          @type theta: positive C{int}          :type sweeps: positive ``int``
716          """          """
717          sweeps=int(sweeps)          sweeps=int(sweeps)
718          if sweeps<1:          if sweeps<1:
# Line 634  class SolverOptions(object): Line 722  class SolverOptions(object):
722          """          """
723          Returns he number of sweeps in the pre-smoothing step of a multi level solver or preconditioner          Returns he number of sweeps in the pre-smoothing step of a multi level solver or preconditioner
724    
725          @rtype: C{int}          :rtype: ``int``
726          """          """
727          return self.__pre_sweeps          return self.__pre_sweeps
728      def setNumPostSweeps(self,sweeps=2):      def setNumPostSweeps(self,sweeps=2):
729          """          """
730          Sets the number of sweeps in the post-smoothing step of a multi level solver or preconditioner          Sets the number of sweeps in the post-smoothing step of a multi level solver or preconditioner
731    
732          @param sweeps: number of sweeps          :param sweeps: number of sweeps
733          @type theta: positive C{int}          :type sweeps: positive ``int``
734          """          """
735          sweeps=int(sweeps)          sweeps=int(sweeps)
736          if sweeps<1:          if sweeps<1:
# Line 652  class SolverOptions(object): Line 740  class SolverOptions(object):
740          """          """
741          Returns he number of sweeps in the post-smoothing step of a multi level solver or preconditioner          Returns he number of sweeps in the post-smoothing step of a multi level solver or preconditioner
742    
743          @rtype: C{int}          :rtype: ``int``
744          """          """
745          return self.__post_sweeps          return self.__post_sweeps
746    
# Line 660  class SolverOptions(object): Line 748  class SolverOptions(object):
748          """          """
749          Sets the relative tolerance for the solver          Sets the relative tolerance for the solver
750    
751          @param rtol: relative tolerance          :param rtol: relative tolerance
752          @type rtol: non-negative C{float}          :type rtol: non-negative ``float``
753          """          """
754          rtol=float(rtol)          rtol=float(rtol)
755          if rtol<0 or rtol>1:          if rtol<0 or rtol>1:
# Line 671  class SolverOptions(object): Line 759  class SolverOptions(object):
759          """          """
760          Returns the relative tolerance for the solver          Returns the relative tolerance for the solver
761    
762          @rtype: C{float}          :rtype: ``float``
763          """          """
764          return self.__tolerance          return self.__tolerance
765      def setAbsoluteTolerance(self,atol=0.):      def setAbsoluteTolerance(self,atol=0.):
766          """          """
767          Sets the absolute tolerance for the solver          Sets the absolute tolerance for the solver
768    
769          @param atol:  absolute tolerance          :param atol:  absolute tolerance
770          @type atol: non-negative C{float}          :type atol: non-negative ``float``
771          """          """
772          atol=float(atol)          atol=float(atol)
773          if atol<0:          if atol<0:
# Line 689  class SolverOptions(object): Line 777  class SolverOptions(object):
777          """          """
778          Returns the absolute tolerance for the solver          Returns the absolute tolerance for the solver
779    
780          @rtype: C{float}          :rtype: ``float``
781          """          """
782          return self.__absolute_tolerance          return self.__absolute_tolerance
783    
784      def setInnerTolerance(self,rtol=0.9):      def setInnerTolerance(self,rtol=0.9):
785          """          """
786           Sets the relative tolerance for an inner iteration scheme for instance           Sets the relative tolerance for an inner iteration scheme for instance on the coarsest level in a multi-level scheme.
         on the coarsest level in a multi-level scheme.  
787    
788          @param rtol: inner relative tolerance          :param rtol: inner relative tolerance
789          @type rtol: positive C{float}          :type rtol: positive ``float``
790          """          """
791          rtol=float(rtol)          rtol=float(rtol)
792          if rtol<=0 or rtol>1:          if rtol<=0 or rtol>1:
# Line 709  class SolverOptions(object): Line 796  class SolverOptions(object):
796          """          """
797          Returns the relative tolerance for an inner iteration scheme          Returns the relative tolerance for an inner iteration scheme
798    
799          @rtype: C{float}          :rtype: ``float``
800          """          """
801          return self.__inner_tolerance          return self.__inner_tolerance
802      def setDropTolerance(self,drop_tol=0.01):      def setDropTolerance(self,drop_tol=0.01):
803          """          """
804          Sets the relative drop tolerance in ILUT          Sets the relative drop tolerance in ILUT
805    
806          @param drop_tol: drop tolerance          :param drop_tol: drop tolerance
807          @type drop_tol: positive C{float}          :type drop_tol: positive ``float``
808          """          """
809          drop_tol=float(drop_tol)          drop_tol=float(drop_tol)
810          if drop_tol<=0 or drop_tol>1:          if drop_tol<=0 or drop_tol>1:
# Line 727  class SolverOptions(object): Line 814  class SolverOptions(object):
814          """          """
815          Returns the relative drop tolerance in ILUT          Returns the relative drop tolerance in ILUT
816    
817          @rtype: C{float}          :rtype: ``float``
818          """          """
819          return self.__drop_tolerance          return self.__drop_tolerance
820      def setDropStorage(self,storage=2.):      def setDropStorage(self,storage=2.):
821          """          """
822          Sets the maximum allowed increase in storage for ILUT. C{storage}=2 would mean that          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.          a doubling of the storage needed for the coefficient matrix is allowed in the ILUT factorization.
824    
825          @param storage: allowed storage increase          :param storage: allowed storage increase
826          @type storage: C{float}          :type storage: ``float``
827          """          """
828          storage=float(storage)          storage=float(storage)
829          if storage<1:          if storage<1:
# Line 747  class SolverOptions(object): Line 834  class SolverOptions(object):
834          """          """
835          Returns the maximum allowed increase in storage for ILUT          Returns the maximum allowed increase in storage for ILUT
836    
837          @rtype: C{float}          :rtype: ``float``
838          """          """
839          return self.__drop_storage          return self.__drop_storage
840      def setRelaxationFactor(self,factor=0.3):      def setRelaxationFactor(self,factor=0.3):
841          """          """
842          Sets the relaxation factor used to add dropped elements in RILU to the main diagonal.          Sets the relaxation factor used to add dropped elements in RILU to the main diagonal.
843    
844          @param factor: relaxation factor          :param factor: relaxation factor
845          @type factor: C{float}          :type factor: ``float``
846          @note: RILU with a relaxation factor 0 is identical to ILU0          :note: RILU with a relaxation factor 0 is identical to ILU0
847          """          """
848          factor=float(factor)          factor=float(factor)
849          if factor<0:          if factor<0:
# Line 767  class SolverOptions(object): Line 854  class SolverOptions(object):
854          """          """
855          Returns the relaxation factor used to add dropped elements in RILU to the main diagonal.          Returns the relaxation factor used to add dropped elements in RILU to the main diagonal.
856    
857          @rtype: C{float}          :rtype: ``float``
858          """          """
859          return self.__relaxation          return self.__relaxation
860      def isSymmetric(self):      def isSymmetric(self):
861          """          """
862          Checks if symmetry of the  coefficient matrix is indicated.          Checks if symmetry of the  coefficient matrix is indicated.
863    
864          @return: True if a symmetric PDE is indicated, False otherwise          :return: True if a symmetric PDE is indicated, False otherwise
865          @rtype: C{bool}          :rtype: ``bool``
866          """          """
867          return self.__symmetric          return self.__symmetric
868      def setSymmetryOn(self):      def setSymmetryOn(self):
# Line 790  class SolverOptions(object): Line 877  class SolverOptions(object):
877          self.__symmetric=False          self.__symmetric=False
878      def setSymmetry(self,flag=False):      def setSymmetry(self,flag=False):
879          """          """
880          Sets the symmetry flag for the coefficient matrix to C{flag}.          Sets the symmetry flag for the coefficient matrix to ``flag``.
881    
882          @param flag: If True, the symmetry flag is set otherwise reset.          :param flag: If True, the symmetry flag is set otherwise reset.
883          @type flag: C{bool}          :type flag: ``bool``
884          """          """
885          if flag:          if flag:
886              self.setSymmetryOn()              self.setSymmetryOn()
# Line 801  class SolverOptions(object): Line 888  class SolverOptions(object):
888              self.setSymmetryOff()              self.setSymmetryOff()
889      def isVerbose(self):      def isVerbose(self):
890          """          """
891          Returns C{True} if the solver is expected to be verbose.          Returns ``True`` if the solver is expected to be verbose.
892    
893          @return: True if verbosity of switched on.          :return: True if verbosity of switched on.
894          @rtype: C{bool}          :rtype: ``bool``
895          """          """
896          return self.__verbose          return self.__verbose
897    
# Line 820  class SolverOptions(object): Line 907  class SolverOptions(object):
907          self.__verbose=False          self.__verbose=False
908      def setVerbosity(self,verbose=False):      def setVerbosity(self,verbose=False):
909          """          """
910          Sets the verbosity flag for the solver to C{flag}.          Sets the verbosity flag for the solver to ``flag``.
911    
912          @param flag: If C{True}, the verbosity of the solver is switched on.          :param verbose: If ``True``, the verbosity of the solver is switched on.
913          @type flag: C{bool}          :type verbose: ``bool``
914          """          """
915          if verbose:          if verbose:
916              self.setVerbosityOn()              self.setVerbosityOn()
# Line 832  class SolverOptions(object): Line 919  class SolverOptions(object):
919                    
920      def adaptInnerTolerance(self):      def adaptInnerTolerance(self):
921          """          """
922          Returns C{True} if the tolerance of the inner solver is selected automatically.          Returns ``True`` if the tolerance of the inner solver is selected automatically.
923          Otherwise the inner tolerance set by L{setInnerTolerance} is used.          Otherwise the inner tolerance set by `setInnerTolerance` is used.
924    
925          @return: C{True} if inner tolerance adaption is chosen.          :return: ``True`` if inner tolerance adaption is chosen.
926          @rtype: C{bool}          :rtype: ``bool``
927          """          """
928          return self.__adapt_inner_tolerance          return self.__adapt_inner_tolerance
929    
# Line 852  class SolverOptions(object): Line 939  class SolverOptions(object):
939          self.__adapt_inner_tolerance=False          self.__adapt_inner_tolerance=False
940      def setInnerToleranceAdaption(self,adapt=True):      def setInnerToleranceAdaption(self,adapt=True):
941          """          """
942          Sets a flag to indicate automatic selection of the inner tolerance.          Sets the flag to indicate automatic selection of the inner tolerance.
943    
944          @param adapt: If C{True}, the inner tolerance is selected automatically.          :param adapt: If ``True``, the inner tolerance is selected automatically.
945          @type adapt: C{bool}          :type adapt: ``bool``
946          """          """
947          if adapt:          if adapt:
948              self.setInnerToleranceAdaptionOn()              self.setInnerToleranceAdaptionOn()
# Line 864  class SolverOptions(object): Line 951  class SolverOptions(object):
951    
952      def acceptConvergenceFailure(self):      def acceptConvergenceFailure(self):
953          """          """
954          Returns C{True} if a failure to meet the stopping criteria within the          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          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          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          continue even if the returned the solution does not necessarily meet the
958          stopping criteria. One can use the L{hasConverged} method to check if the          stopping criteria. One can use the `hasConverged` method to check if the
959          last call to the solver was successful.          last call to the solver was successful.
960    
961          @return: C{True} if a failure to achieve convergence is accepted.          :return: ``True`` if a failure to achieve convergence is accepted.
962          @rtype: C{bool}          :rtype: ``bool``
963          """          """
964          return self.__accept_convergence_failure          return self.__accept_convergence_failure
965    
# Line 888  class SolverOptions(object): Line 975  class SolverOptions(object):
975          self.__accept_convergence_failure=False          self.__accept_convergence_failure=False
976      def setAcceptanceConvergenceFailure(self,accept=False):      def setAcceptanceConvergenceFailure(self,accept=False):
977          """          """
978          Sets a flag to indicate the acceptance of a failure of convergence.          Sets the flag to indicate the acceptance of a failure of convergence.
979    
980          @param accept: If C{True}, any failure to achieve convergence is accepted.          :param accept: If ``True``, any failure to achieve convergence is accepted.
981          @type accept: C{bool}          :type accept: ``bool``
982          """          """
983          if accept:          if accept:
984              self.setAcceptanceConvergenceFailureOn()              self.setAcceptanceConvergenceFailureOn()
985          else:          else:
986              self.setAcceptanceConvergenceFailureOff()              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 927  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 979  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 1012  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 1049  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 1058  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 1099  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 1112  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 1126  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 1174  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 1186  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 1199  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 1213  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 1238  class PDECoef(object): Line 1509  class PDECoef(object):
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*
1520    
1521     where M{L} is an operator and M{f} is the right hand side. This operator     where *L* is an operator and *f* is the right hand side. This operator
1522     problem will be solved to get the unknown M{u}.     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 1292  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 1313  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 1327  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 1339  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 1356  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):     def getDomainStatus(self):
# Line 1373  class LinearProblem(object): Line 1644  class LinearProblem(object):
1644       return self.__system_status       return self.__system_status
1645     def setSystemStatus(self,status=None):     def setSystemStatus(self,status=None):
1646       """       """
1647       Sets the system status to C{status} if C{status} is not present the       Sets the system status to ``status`` if ``status`` is not present the
1648       current status of the domain is used.       current status of the domain is used.
1649       """       """
1650       if status == None:       if status == None:
# Line 1385  class LinearProblem(object): Line 1656  class LinearProblem(object):
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 1394  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 1409  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 1424  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 1434  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 1455  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 1473  class LinearProblem(object): Line 1744  class LinearProblem(object):
1744         """         """
1745         Sets the solver options.         Sets the solver options.
1746    
1747         @param options: the new solver options. If equal C{None}, the solver options are set to the default.         :param options: the new solver options. If equal ``None``, the solver options are set to the default.
1748         @type options: L{SolverOptions} or C{None}         :type options: `SolverOptions` or ``None``
1749         @note: The symmetry flag of options is overwritten by the symmetry flag of the L{LinearProblem}.         :note: The symmetry flag of options is overwritten by the symmetry flag of the `LinearProblem`.
1750         """         """
1751         if options==None:         if options==None:
1752            self.__solver_options=SolverOptions()            self.__solver_options=SolverOptions()
# Line 1489  class LinearProblem(object): Line 1760  class LinearProblem(object):
1760         """         """
1761         Returns the solver options         Returns the solver options
1762        
1763         @rtype: L{SolverOptions}         :rtype: `SolverOptions`
1764         """         """
1765         self.__solver_options.setSymmetry(self.__sym)         self.__solver_options.setSymmetry(self.__sym)
1766         return self.__solver_options         return self.__solver_options
# Line 1498  class LinearProblem(object): Line 1769  class LinearProblem(object):
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.getSolverOptions().getSolverMethod()==self.getSolverOptions().LUMPING        return self.getSolverOptions().getSolverMethod() in [ SolverOptions.ROWSUM_LUMPING, SolverOptions.HRZ_LUMPING ]
1776     # ==========================================================================     # ==========================================================================
1777     #    symmetry  flag:     #    symmetry  flag:
1778     # ==========================================================================     # ==========================================================================
# Line 1509  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()        :note: the method is equivalent to use getSolverOptions().isSymmetric()
1786        """        """
1787        self.getSolverOptions().isSymmetric()        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        :note: The method overwrites the symmetry flag set by the solver options
1793        """        """
1794        self.__sym=True        self.__sym=True
1795        self.getSolverOptions().setSymmetryOn()        self.getSolverOptions().setSymmetryOn()
# Line 1526  class LinearProblem(object): Line 1797  class LinearProblem(object):
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        :note: The method overwrites the symmetry flag set by the solver options
1801        """        """
1802        self.__sym=False        self.__sym=False
1803        self.getSolverOptions().setSymmetryOff()        self.getSolverOptions().setSymmetryOff()
1804    
1805     def setSymmetry(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        :note: The method overwrites the symmetry flag set by the solver options
1812        """        """
1813        self.getSolverOptions().setSymmetry(flag)        self.getSolverOptions().setSymmetry(flag)
1814     # ==========================================================================     # ==========================================================================
# Line 1547  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 1557  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 1568  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 1583  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 1597  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 1611  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 1627  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 1641  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 1655  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:
# Line 1676  class LinearProblem(object): Line 1947  class LinearProblem(object):
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.        SMALL_TOLERANCE=util.EPSILON*10.
1959        A=self.getCoefficient(name)        A=self.getCoefficient(name)
# Line 1723  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.        SMALL_TOLERANCE=util.EPSILON*10.
2009        B=self.getCoefficient(name0)        B=self.getCoefficient(name0)
# Line 1781  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 1796  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 1822  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 1838  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 1860  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 1972  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 1980  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 2015  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 with the tolrance         Sets the solution assuming that makes the system valid with the tolrance
2292         defined by the solver options         defined by the solver options
2293         """         """
2294         self.__solution_rtol=self.getSolverOptions().getTolerance()         if validate:
2295         self.__solution_atol=self.getSolverOptions().getAbsoluteTolerance()        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 2064  class LinearProblem(object): Line 2335  class LinearProblem(object):
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):
# Line 2077  class LinearProblem(object): Line 2352  class LinearProblem(object):
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 2136  class LinearProblem(object): Line 2411  class LinearProblem(object):
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 None        return None
# Line 2145  class LinearProblem(object): Line 2420  class LinearProblem(object):
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 2154  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 2169  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 2181  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())
# Line 2192  class LinearProblem(object): Line 2467  class LinearProblem(object):
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 2325  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 2367  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    
# Line 2376  class LinearPDE(LinearProblem): Line 2651  class LinearPDE(LinearProblem):
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        solver_options=self.getSolverOptions()        if self.isUsingLumping():
2655        return self.getDomain().getSystemMatrixTypeId(solver_options.getSolverMethod(), solver_options.getPreconditioner(),solver_options.getPackage(), solver_options.isSymmetric())       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 2421  class LinearPDE(LinearProblem): Line 2699  class LinearPDE(LinearProblem):
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`
2704         """         """
2705         option_class=self.getSolverOptions()         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:
2713               self.trace("PDE is resolved.")               self.trace("PDE is resolved.")
# Line 2439  class LinearPDE(LinearProblem): Line 2719  class LinearPDE(LinearProblem):
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 2500  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 2616  class LinearPDE(LinearProblem): Line 2897  class LinearPDE(LinearProblem):
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 2646  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 2735  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 2749  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 2775  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 2793  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 2808  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 2846  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 2864  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 2884  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():
# Line 2948  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 2964  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 2982  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 3053  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 3068  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 3080  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 3219  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 3272  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 3289  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 3318  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``
3608        @type M_reduced: any type that can be cast to a L{Data<escript.Data>}        :type M_reduced: any type that can be cast to a `Data`
3609                         object on L{Function<escript.ReducedFunction>}                         object on `Function`
3610        @keyword A: value for coefficient C{A}        :keyword A: value for coefficient ``A``
3611        @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
3612                 L{Function<escript.Function>}                 `Function`
3613        @keyword A_reduced: value for coefficient C{A_reduced}        :keyword A_reduced: value for coefficient ``A_reduced``
3614        @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`
3615                         object on L{ReducedFunction<escript.ReducedFunction>}                         object on `ReducedFunction`
3616        @keyword B: value for coefficient C{B}        :keyword B: value for coefficient ``B``
3617        @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
3618                 L{Function<escript.Function>}                 `Function`
3619        @keyword B_reduced: value for coefficient C{B_reduced}        :keyword B_reduced: value for coefficient ``B_reduced``
3620        @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`
3621                         object on L{ReducedFunction<escript.ReducedFunction>}                         object on `ReducedFunction`
3622        @keyword C: value for coefficient C{C}        :keyword C: value for coefficient ``C``
3623        @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
3624                 L{Function<escript.Function>}                 `Function`
3625        @keyword C_reduced: value for coefficient C{C_reduced}        :keyword C_reduced: value for coefficient ``C_reduced``
3626        @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`
3627                         object on L{ReducedFunction<escript.ReducedFunction>}                         object on `ReducedFunction`
3628        @keyword D: value for coefficient C{D}        :keyword D: value for coefficient ``D``
3629        @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
3630                 L{Function<escript.Function>}                 `Function`
3631        @keyword D_reduced: value for coefficient C{D_reduced}        :keyword D_reduced: value for coefficient ``D_reduced``
3632        @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`
3633                         object on L{ReducedFunction<escript.ReducedFunction>}                         object on `ReducedFunction`
3634        @keyword X: value for coefficient C{X}        :keyword X: value for coefficient ``X``
3635        @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
3636                 L{Function<escript.Function>}                 `Function`
3637        @keyword X_reduced: value for coefficient C{X_reduced}        :keyword X_reduced: value for coefficient ``X_reduced``
3638        @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`
3639                         object on L{ReducedFunction<escript.ReducedFunction>}                         object on `ReducedFunction`
3640        @keyword Y: value for coefficient C{Y}        :keyword Y: value for coefficient ``Y``
3641        @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
3642                 L{Function<escript.Function>}                 `Function`
3643        @keyword Y_reduced: value for coefficient C{Y_reduced}        :keyword Y_reduced: value for coefficient ``Y_reduced``
3644        @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`
3645                         object on L{ReducedFunction<escript.Function>}                         object on `ReducedFunction`
3646        @keyword m: value for coefficient C{m}        :keyword m: value for coefficient ``m``
3647        @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
3648                 L{FunctionOnBoundary<escript.FunctionOnBoundary>}                 `FunctionOnBoundary`
3649        @keyword m_reduced: value for coefficient C{m_reduced}        :keyword m_reduced: value for coefficient ``m_reduced``
3650        @type m_reduced: any type that can be cast to a L{Data<escript.Data>}        :type m_reduced: any type that can be cast to a `Data`
3651                         object on L{FunctionOnBoundary<escript.ReducedFunctionOnBoundary>}                         object on `FunctionOnBoundary`
3652        @keyword d: value for coefficient C{d}        :keyword d: value for coefficient ``d``
3653        @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
3654                 L{FunctionOnBoundary<escript.FunctionOnBoundary>}                 `FunctionOnBoundary`
3655        @keyword d_reduced: value for coefficient C{d_reduced}        :keyword d_reduced: value for coefficient ``d_reduced``
3656        @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`
3657                         object on L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}                         object on `ReducedFunctionOnBoundary`
3658        @keyword y: value for coefficient C{y}        :keyword y: value for coefficient ``y``
3659        @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
3660                 L{FunctionOnBoundary<escript.FunctionOnBoundary>}                 `FunctionOnBoundary`
3661        @keyword d_contact: value for coefficient C{d_contact}        :keyword d_contact: value for coefficient ``d_contact``
3662        @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`
3663                         object on L{FunctionOnContactOne<escript.FunctionOnContactOne>} or L{FunctionOnContactZero<escript.FunctionOnContactZero>}                         object on `FunctionOnContactOne` or `FunctionOnContactZero`
3664        @keyword d_contact_reduced: value for coefficient C{d_contact_reduced}        :keyword d_contact_reduced: value for coefficient ``d_contact_reduced``
3665        @type d_contact_reduced: any type that can be cast to a L{Data<escript.Data>} object on L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>} or L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>}        :type d_contact_reduced: any type that can be cast to a `Data` object on `ReducedFunctionOnContactOne` or `ReducedFunctionOnContactZero`
3666        @keyword y_contact: value for coefficient C{y_contact}        :keyword y_contact: value for coefficient ``y_contact``
3667        @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`
3668                         object on L{FunctionOnContactOne<escript.FunctionOnContactOne>} or L{FunctionOnContactZero<escript.FunctionOnContactZero>}                         object on `FunctionOnContactOne` or `FunctionOnContactZero`
3669        @keyword y_contact_reduced: value for coefficient C{y_contact_reduced}        :keyword y_contact_reduced: value for coefficient ``y_contact_reduced``
3670        @type y_contact_reduced: any type that can be cast to a L{Data<escript.Data>} object on L{ReducedFunctionOnContactOne<escript.FunctionOnContactOne>} or L{ReducedFunctionOnContactZero<escript.FunctionOnContactZero>}        :type y_contact_reduced: any type that can be cast to a `Data` object on `ReducedFunctionOnContactOne` or `ReducedFunctionOnContactZero`
3671        @keyword r: values prescribed to the solution at the locations of constraints        :keyword r: values prescribed to the solution at the locations of constraints
3672        @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
3673                 L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}                 `Solution` or `ReducedSolution`
3674                 depending on whether reduced order is used for the solution                 depending on whether reduced order is used for the solution
3675        @keyword q: mask for the location of constraints        :keyword q: mask for the location of constraints
3676        @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
3677                 L{Solution<escript.Solution>} or                 `Solution` or
3678                 L{ReducedSolution<escript.ReducedSolution>} depending on whether                 `ReducedSolution` depending on whether
3679                 reduced order is used for the representation of the equation                 reduced order is used for the representation of the equation
3680        @raise IllegalCoefficient: if an unknown coefficient keyword is used        :raise IllegalCoefficient: if an unknown coefficient keyword is used
3681        """        """
3682        super(TransportPDE,self).setValue(**coefficients)        super(TransportPDE,self).setValue(**coefficients)
3683    
# Line 3403  class TransportPDE(LinearProblem): Line 3685  class TransportPDE(LinearProblem):
3685         """         """
3686         Returns an instance of a new transport operator.         Returns an instance of a new transport operator.
3687         """         """
        if self.useBackwardEuler():  
          theta=1.  
        else:  
          theta=0.5  
3688         optype=self.getRequiredOperatorType()         optype=self.getRequiredOperatorType()
3689         self.trace("New Transport problem pf type %s is allocated."%optype)         self.trace("New Transport problem pf type %s is allocated."%optype)
3690         return self.getDomain().newTransportProblem( \         return self.getDomain().newTransportProblem( \
3691                                 theta,                                 self.useBackwardEuler(),
3692                                 self.getNumEquations(), \                                 self.getNumEquations(), \
3693                                 self.getFunctionSpaceForSolution(), \                                 self.getFunctionSpaceForSolution(), \
3694                                 optype)                                 optype)
3695    
    def setInitialSolution(self,u):  
        """  
        Sets the initial solution.  
   
        @param u: new initial solution  
        @type u: any object that can be interpolated to a L{Data<escript.Data>}  
                 object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}  
        @note: C{u} must be non-negative  
        """  
        u2=util.interpolate(u,self.getFunctionSpaceForSolution())  
        if self.getNumSolutions() == 1:  
           if u2.getShape()!=():  
               raise ValueError,"Illegal shape %s of initial solution."%(u2.getShape(),)  
        else:  
           if u2.getShape()!=(self.getNumSolutions(),):  
               raise ValueError,"Illegal shape %s of initial solution."%(u2.getShape(),)  
        self.getOperator().setInitialValue(u2)  
3696    
3697     def getRequiredOperatorType(self):     def getRequiredOperatorType(self):
3698        """        """
3699        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.
3700    
3701        @return: a code to indicate the type of transport problem scheme used        :return: a code to indicate the type of transport problem scheme used
3702        @rtype: C{float}        :rtype: ``float``
3703        """        """
3704        solver_options=self.getSolverOptions()        solver_options=self.getSolverOptions()
3705        return self.getDomain().getTransportTypeId(solver_options.getSolverMethod(), solver_options.getPreconditioner(),solver_options.getPackage(), solver_options.isSymmetric())        return self.getDomain().getTransportTypeId(solver_options.getSolverMethod(), solver_options.getPreconditioner(),solver_options.getPackage(), solver_options.isSymmetric())
3706    
3707     def getUnlimitedTimeStepSize(self):     def getUnlimitedTimeStepSize(self):
3708        """        """
3709        Returns the value returned by the C{getSafeTimeStepSize} method to        Returns the value returned by the ``getSafeTimeStepSize`` method to
3710        indicate no limit on the safe time step size.        indicate no limit on the safe time step size.
3711    
3712         @return: the value used to indicate that no limit is set to the time         :return: the value used to indicate that no limit is set to the time
3713                  step size                  step size
3714         @rtype: C{float}         :rtype: ``float``
3715         @note: Typically the biggest positive float is returned         :note: Typically the biggest positive float is returned
3716        """        """
3717        return self.getOperator().getUnlimitedTimeStepSize()        return self.getOperator().getUnlimitedTimeStepSize()
3718    
# Line 3459  class TransportPDE(LinearProblem): Line 3720  class TransportPDE(LinearProblem):
3720         """         """
3721         Returns a safe time step size to do the next time step.         Returns a safe time step size to do the next time step.
3722    
3723         @return: safe time step size         :return: safe time step size
3724         @rtype: C{float}         :rtype: ``float``
3725         @note: If not C{getSafeTimeStepSize()} < C{getUnlimitedTimeStepSize()}         :note: If not ``getSafeTimeStepSize()`` < ``getUnlimitedTimeStepSize()``
3726                any time step size can be used.                any time step size can be used.
3727         """         """
3728         return self.getOperator().getSafeTimeStepSize()         return self.getOperator().getSafeTimeStepSize()
# Line 3470  class TransportPDE(LinearProblem): Line 3731  class TransportPDE(LinearProblem):
3731         """         """
3732         Sets the weighting factor used to insert the constraints into the problem         Sets the weighting factor used to insert the constraints into the problem
3733    
3734         @param value: value for the weighting factor         :param value: value for the weighting factor
3735         @type value: large positive C{float}         :type value: large positive ``float``
3736         """         """
3737         if not value>0:         if not value>0:
3738           raise ValueError,"weighting factor needs to be positive."           raise ValueError,"weighting factor needs to be positive."
# Line 3481  class TransportPDE(LinearProblem): Line 3742  class TransportPDE(LinearProblem):
3742     def getConstraintWeightingFactor(self):     def getConstraintWeightingFactor(self):
3743         """         """
3744         returns the weighting factor used to insert the constraints into the problem         returns the weighting factor used to insert the constraints into the problem
3745         @return: value for the weighting factor         :return: value for the weighting factor
3746         @rtype: C{float}         :rtype: ``float``
3747         """         """
3748         return self.__constraint_factor         return self.__constraint_factor
3749     #====================================================================     #====================================================================
3750     def getSolution(self,dt):     def getSolution(self, dt=None, u0=None):
3751         """         """
3752         Returns the solution of the problem.         Returns the solution by marching forward by time step dt. if ''u0'' is present,
3753           ''u0'' is used as the initial value otherwise the solution from the last call is used.
3754    
3755         @return: the solution         :param dt: time step size. If ``None`` the last solution is returned.
3756         @rtype: L{Data<escript.Data>}         :type dt: positive ``float`` or ``None``
3757         """         :param u0: new initial solution or ``None``
3758         option_class=self.getSolverOptions()         :type u0: any object that can be interpolated to a `Data`
3759         if dt<=0:                  object on `Solution` or `ReducedSolution`
3760             raise ValueError,"step size needs to be positive."         :return: the solution
3761         self.setSolution(self.getOperator().solve(self.getRightHandSide(),dt,option_class))         :rtype: `Data`
3762         self.validSolution()         """
3763           if not dt == None:
3764          option_class=self.getSolverOptions()
3765          if dt<=0:
3766              raise ValueError,"step size needs to be positive."
3767          if u0 == None:
3768              u0=self.getCurrentSolution()
3769          else:
3770              u0=util.interpolate(u0,self.getFunctionSpaceForSolution())
3771              if self.getNumSolutions() == 1:
3772            if u0.getShape()!=():
3773                raise ValueError,"Illegal shape %s of initial solution."%(u0.getShape(),)
3774              else:
3775                if u0.getShape()!=(self.getNumSolutions(),):
3776                  raise ValueError,"Illegal shape %s of initial solution."%(u0.getShape(),)
3777          self.setSolution(self.getOperator().solve(u0, self.getRightHandSide(),dt,option_class))
3778          self.validSolution()
3779         return self.getCurrentSolution()         return self.getCurrentSolution()
3780    
3781       def setInitialSolution(self,u):
3782           """
3783           Sets the initial solution.
3784    
3785           :param u: initial solution
3786           :type u: any object that can be interpolated to a `Data`
3787                    object on `Solution` or `ReducedSolution`
3788           """
3789           u2=util.interpolate(u,self.getFunctionSpaceForSolution())
3790           if self.getNumSolutions() == 1:
3791              if u2.getShape()!=():
3792                  raise ValueError,"Illegal shape %s of initial solution."%(u2.getShape(),)
3793           else:
3794              if u2.getShape()!=(self.getNumSolutions(),):
3795                  raise ValueError,"Illegal shape %s of initial solution."%(u2.getShape(),)
3796           self.setSolution(u2,validate=False)
3797    
3798    
3799     def getSystem(self):     def getSystem(self):
3800         """         """
3801         Returns the operator and right hand side of the PDE.         Returns the operator and right hand side of the PDE.
3802    
3803         @return: the discrete version of the PDE         :return: the discrete version of the PDE
3804         @rtype: C{tuple} of L{Operator,<escript.Operator>} and         :rtype: ``tuple`` of `Operator` and
3805                 L{Data<escript.Data>}                 `Data`
3806    
3807         """         """
3808         if not self.isOperatorValid() or not self.isRightHandSideValid():         if not self.isOperatorValid() or not self.isRightHandSideValid():
# Line 3552  class TransportPDE(LinearProblem): Line 3848  class TransportPDE(LinearProblem):
3848    
3849     def setDebug(self, flag):     def setDebug(self, flag):
3850       """       """
3851       Switches debug output on if C{flag} is True,       Switches debug output on if ``flag`` is True,
3852       otherwise it is switched off.       otherwise it is switched off.
3853    
3854       @param flag: desired debug status       :param flag: desired debug status
3855       @type flag: C{bool}       :type flag: ``bool``
3856       """       """
3857       if flag:       if flag:
3858           self.setDebugOn()           self.setDebugOn()
# Line 3575  class TransportPDE(LinearProblem): Line 3871  class TransportPDE(LinearProblem):
3871       """       """
3872       super(TransportPDE,self).setDebugOff()       super(TransportPDE,self).setDebugOff()
3873            
3874    def SingleTransportPDE(domain,useBackwardEuler=False, debug=False):
3875       """
3876       Defines a single transport problem
3877    
3878       :param domain: domain of the PDE
3879       :type domain: `Domain`
3880       :param debug: if True debug information is printed
3881       :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.
3882       :rtype: `TransportPDE`
3883       """
3884       return TransportPDE(domain,numEquations=1,numSolutions=1, useBackwardEuler=useBackwardEuler, debug=debug)

Legend:
Removed from v.2548  
changed lines
  Added in v.3402

  ViewVC Help
Powered by ViewVC 1.1.26