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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2596 - (hide annotations)
Thu Aug 6 03:09:03 2009 UTC (10 years ago) by lgao
File MIME type: text/x-python
File size: 154269 byte(s)
Add new diagnostic information: "net time" (the execution time of solver, excluding the setup time for solver and the execution time of preconditioner)


1 ksteube 1809
2     ########################################################
3 ksteube 1312 #
4 jfenwick 2548 # Copyright (c) 2003-2009 by University of Queensland
5 ksteube 1809 # Earth Systems Science Computational Center (ESSCC)
6     # http://www.uq.edu.au/esscc
7 ksteube 1312 #
8 ksteube 1809 # Primary Business: Queensland, Australia
9     # Licensed under the Open Software License version 3.0
10     # http://www.opensource.org/licenses/osl-3.0.php
11 ksteube 1312 #
12 ksteube 1809 ########################################################
13 ksteube 1312
14 jfenwick 2549 __copyright__="""Copyright (c) 2003-2009 by University of Queensland
15 ksteube 1809 Earth Systems Science Computational Center (ESSCC)
16     http://www.uq.edu.au/esscc
17     Primary Business: Queensland, Australia"""
18     __license__="""Licensed under the Open Software License version 3.0
19     http://www.opensource.org/licenses/osl-3.0.php"""
20 jfenwick 2344 __url__="https://launchpad.net/escript-finley"
21 ksteube 1809
22 jgs 149 """
23 jgs 151 The module provides an interface to define and solve linear partial
24 caltinay 2169 differential equations (PDEs) and Transport problems within L{escript}.
25     L{linearPDEs} does not provide any solver capabilities in itself but hands the
26     PDE over to the PDE solver library defined through the L{Domain<escript.Domain>}
27     of the PDE. The general interface is provided through the L{LinearPDE} class.
28     L{TransportProblem} provides an interface to initial value problems dominated
29     by its advective terms.
30 jgs 102
31 jgs 149 @var __author__: name of author
32 gross 637 @var __copyright__: copyrights
33 elspeth 614 @var __license__: licence agreement
34 jgs 149 @var __url__: url entry point on documentation
35     @var __version__: version
36     @var __date__: date of the version
37 jgs 102 """
38    
39 gross 1417 import math
40 jgs 102 import escript
41     import util
42 jfenwick 2455 import numpy
43 jgs 102
44 jgs 149 __author__="Lutz Gross, l.gross@uq.edu.au"
45 jgs 102
46 jgs 149
47 gross 2470 class SolverOptions(object):
48     """
49     this class defines the solver options for a linear or non-linear solver.
50    
51     The option also supports the handling of diagnostic informations.
52    
53     Typical usage is
54    
55     opts=SolverOptions()
56     print opts
57     opts.resetDiagnostics()
58     u=solver(opts)
59     print "number of iteration steps: =",opts.getDiagnostics("num_iter")
60    
61    
62     @cvar DEFAULT: The default method used to solve the system of linear equations
63     @cvar DIRECT: The direct solver based on LDU factorization
64     @cvar CHOLEVSKY: The direct solver based on LDLt factorization (can only be applied for symmetric PDEs)
65     @cvar PCG: The preconditioned conjugate gradient method (can only be applied for symmetric PDEs)
66     @cvar CR: The conjugate residual method
67     @cvar CGS: The conjugate gradient square method
68     @cvar BICGSTAB: The stabilized Bi-Conjugate Gradient method
69     @cvar TFQMR: Transport Free Quasi Minimal Residual method
70     @cvar MINRES: Minimum residual method
71     @cvar SSOR: The symmetric over-relaxation method
72     @cvar ILU0: The incomplete LU factorization preconditioner with no fill-in
73     @cvar ILUT: The incomplete LU factorization preconditioner with fill-in
74     @cvar JACOBI: The Jacobi preconditioner
75     @cvar GMRES: The Gram-Schmidt minimum residual method
76     @cvar PRES20: Special GMRES with restart after 20 steps and truncation after 5 residuals
77     @cvar LUMPING: Matrix lumping
78     @cvar NO_REORDERING: No matrix reordering allowed
79     @cvar MINIMUM_FILL_IN: Reorder matrix to reduce fill-in during factorization
80     @cvar NESTED_DISSECTION: Reorder matrix to improve load balancing during factorization
81     @cvar PASO: PASO solver package
82     @cvar SCSL: SGI SCSL solver library
83     @cvar MKL: Intel's MKL solver library
84     @cvar UMFPACK: The UMFPACK library
85     @cvar TRILINOS: The TRILINOS parallel solver class library from Sandia National Labs
86     @cvar ITERATIVE: The default iterative solver
87     @cvar AMG: Algebraic Multi Grid
88     @cvar REC_ILU: recursive ILU0
89     @cvar RILU: relaxed ILU0
90     @cvar GAUSS_SEIDEL: Gauss-Seidel solver
91     @cvar DEFAULT_REORDERING: the reordering method recommended by the solver
92     @cvar SUPER_LU: the Super_LU solver package
93     @cvar PASTIX: the Pastix direct solver_package
94     @cvar YAIR_SHAPIRA_COARSENING: AMG coarsening method by Yair-Shapira
95     @cvar RUGE_STUEBEN_COARSENING: AMG coarsening method by Ruge and Stueben
96     @cvar AGGREGATION_COARSENING: AMG coarsening using (symmetric) aggregation
97 artak 2524 @cvar MIN_COARSE_MATRIX_SIZE: minimum size of the coarsest level matrix to use direct solver.
98 gross 2470 @cvar NO_PRECONDITIONER: no preconditioner is applied.
99     """
100     DEFAULT= 0
101     DIRECT= 1
102     CHOLEVSKY= 2
103     PCG= 3
104     CR= 4
105     CGS= 5
106     BICGSTAB= 6
107     SSOR= 7
108     ILU0= 8
109     ILUT= 9
110     JACOBI= 10
111     GMRES= 11
112     PRES20= 12
113     LUMPING= 13
114     NO_REORDERING= 17
115     MINIMUM_FILL_IN= 18
116     NESTED_DISSECTION= 19
117     MKL= 15
118     UMFPACK= 16
119     ITERATIVE= 20
120     PASO= 21
121     AMG= 22
122     REC_ILU = 23
123     TRILINOS = 24
124     NONLINEAR_GMRES = 25
125     TFQMR = 26
126     MINRES = 27
127     GAUSS_SEIDEL=28
128     RILU=29
129     DEFAULT_REORDERING=30
130     SUPER_LU=31
131     PASTIX=32
132     YAIR_SHAPIRA_COARSENING=33
133     RUGE_STUEBEN_COARSENING=34
134     AGGREGATION_COARSENING=35
135     NO_PRECONDITIONER=36
136 artak 2524 MIN_COARSE_MATRIX_SIZE=37
137    
138 gross 2470 def __init__(self):
139     self.setLevelMax()
140     self.setCoarseningThreshold()
141     self.setNumSweeps()
142     self.setNumPreSweeps()
143     self.setNumPostSweeps()
144     self.setTolerance()
145     self.setAbsoluteTolerance()
146     self.setInnerTolerance()
147     self.setDropTolerance()
148     self.setDropStorage()
149     self.setIterMax()
150     self.setInnerIterMax()
151     self.setTruncation()
152     self.setRestart()
153     self.setSymmetry()
154     self.setVerbosity()
155     self.setInnerToleranceAdaption()
156     self.setAcceptanceConvergenceFailure()
157     self.setReordering()
158     self.setPackage()
159     self.setSolverMethod()
160     self.setPreconditioner()
161     self.setCoarsening()
162 artak 2524 self.setMinCoarseMatrixSize()
163 gross 2470 self.setRelaxationFactor()
164     self.resetDiagnostics(all=True)
165    
166     def __str__(self):
167     return self.getSummary()
168     def getSummary(self):
169     """
170     Returns a string reporting the current settings
171     """
172     out="Solver Package: %s"%(self.getName(self.getPackage()))
173     out+="\nVerbosity = %s"%self.isVerbose()
174     out+="\nAccept failed convergence = %s"%self.acceptConvergenceFailure()
175     out+="\nRelative tolerance = %e"%self.getTolerance()
176     out+="\nAbsolute tolerance = %e"%self.getAbsoluteTolerance()
177     out+="\nSymmetric problem = %s"%self.isSymmetric()
178     out+="\nMaximum number of iteration steps = %s"%self.getIterMax()
179 gross 2474 out+="\nInner tolerance = %e"%self.getInnerTolerance()
180     out+="\nAdapt innner tolerance = %s"%self.adaptInnerTolerance()
181    
182 gross 2470 if self.getPackage() == self.PASO:
183     out+="\nSolver method = %s"%self.getName(self.getSolverMethod())
184     if self.getSolverMethod() == self.GMRES:
185     out+="\nTruncation = %s"%self.getTruncation()
186     out+="\nRestart = %s"%self.getRestart()
187     if self.getSolverMethod() == self.AMG:
188     out+="\nNumber of pre / post sweeps = %s / %s, %s"%(self.getNumPreSweeps(), self.getNumPostSweeps(), self.getNumSweeps())
189     out+="\nMaximum number of levels = %s"%self.LevelMax()
190     out+="\nCoarsening threshold = %e"%self.getCoarseningThreshold()
191     out+="\Coarsening method = %s"%self.getName(self.getCoarsening())
192     out+="\nPreconditioner = %s"%self.getName(self.getPreconditioner())
193     if self.getPreconditioner() == self.AMG:
194     out+="\nMaximum number of levels = %s"%self.LevelMax()
195     out+="\nCoarsening method = %s"%self.getName(self.getCoarsening())
196 artak 2524 out+="\nCoarsening threshold = %e"%self.getMinCoarseMatrixSize()
197     out+="\nMinimum size of the coarsest level matrix = %e"%self.getCoarseningThreshold()
198 gross 2470 out+="\nNumber of pre / post sweeps = %s / %s, %s"%(self.getNumPreSweeps(), self.getNumPostSweeps(), self.getNumSweeps())
199     if self.getPreconditioner() == self.GAUSS_SEIDEL:
200     out+="\nNumber of sweeps = %s"%self.getNumSweeps()
201     if self.getPreconditioner() == self.ILUT:
202     out+="\nDrop tolerance = %e"%self.getDropTolerance()
203     out+="\nStorage increase = %e"%self.getDropStorage()
204     if self.getPreconditioner() == self.RILU:
205     out+="\nRelaxation factor = %e"%self.getRelaxationFactor()
206     return out
207    
208     def getName(self,key):
209     """
210     returns the name of a given key
211    
212     @param key: a valid key
213     """
214     if key == self.DEFAULT: return "DEFAULT"
215     if key == self.DIRECT: return "DIRECT"
216     if key == self.CHOLEVSKY: return "CHOLEVSKY"
217     if key == self.PCG: return "PCG"
218     if key == self.CR: return "CR"
219     if key == self.CGS: return "CGS"
220     if key == self.BICGSTAB: return "BICGSTAB"
221     if key == self.SSOR: return "SSOR"
222     if key == self.ILU0: return "ILU0:"
223     if key == self.ILUT: return "ILUT"
224     if key == self.JACOBI: return "JACOBI"
225     if key == self.GMRES: return "GMRES"
226     if key == self.PRES20: return "PRES20"
227     if key == self.LUMPING: return "LUMPING"
228     if key == self.NO_REORDERING: return "NO_REORDERING"
229     if key == self.MINIMUM_FILL_IN: return "MINIMUM_FILL_IN"
230     if key == self.NESTED_DISSECTION: return "NESTED_DISSECTION"
231     if key == self.MKL: return "MKL"
232     if key == self.UMFPACK: return "UMFPACK"
233     if key == self.ITERATIVE: return "ITERATIVE"
234     if key == self.PASO: return "PASO"
235     if key == self.AMG: return "AMG"
236     if key == self.REC_ILU: return "REC_ILU"
237     if key == self.TRILINOS: return "TRILINOS"
238     if key == self.NONLINEAR_GMRES: return "NONLINEAR_GMRES"
239     if key == self.TFQMR: return "TFQMR"
240     if key == self.MINRES: return "MINRES"
241     if key == self.GAUSS_SEIDEL: return "GAUSS_SEIDEL"
242     if key == self.RILU: return "RILU"
243     if key == self.DEFAULT_REORDERING: return "DEFAULT_REORDERING"
244     if key == self.SUPER_LU: return "SUPER_LU"
245     if key == self.PASTIX: return "PASTIX"
246     if key == self.YAIR_SHAPIRA_COARSENING: return "YAIR_SHAPIRA_COARSENING"
247     if key == self.RUGE_STUEBEN_COARSENING: return "RUGE_STUEBEN_COARSENING"
248     if key == self.AGGREGATION_COARSENING: return "AGGREGATION_COARSENING"
249     if key == self.NO_PRECONDITIONER: return "NO_PRECONDITIONER"
250 artak 2524 if key == self.MIN_COARSE_MATRIX_SIZE: return "MIN_COARSE_MATRIX_SIZE"
251 gross 2470
252     def resetDiagnostics(self,all=False):
253     """
254     resets the diagnostics
255    
256     @param all: if C{all} is C{True} all diagnostics including accumulative counters are reset.
257     @type all: C{bool}
258     """
259     self.__num_iter=None
260     self.__num_level=None
261     self.__num_inner_iter=None
262     self.__time=None
263     self.__set_up_time=None
264 lgao 2596 self.__net_time=None
265 gross 2470 self.__residual_norm=None
266     self.__converged=None
267     if all:
268     self.__cum_num_inner_iter=0
269     self.__cum_num_iter=0
270     self.__cum_time=0
271     self.__cum_set_up_time=0
272 lgao 2596 self.__cum_net_time=0
273 gross 2470
274     def _updateDiagnostics(self, name, value):
275     """
276     Updates diagnostic information
277    
278     @param name: name of diagnostic information
279 lgao 2596 @type name: C{str} in the list "num_iter", "num_level", "num_inner_iter", "time", "set_up_time", "net_time", "residual_norm", "converged".
280 gross 2470 @param vale: new value of the diagnostic information
281     @note: this function is used by a solver to report diagnostics informations.
282     """
283     if name == "num_iter":
284     self.__num_iter=int(value)
285     self.__cum_num_iter+=self.__num_iter
286     if name == "num_level":
287 lgao 2596 self.__num_level=int(value)
288 gross 2470 if name == "num_inner_iter":
289     self.__num_inner_iter=int(value)
290     self.__cum_num_inner_iter+=self.__num_inner_iter
291     if name == "time":
292     self.__time=float(value)
293     self.__cum_time+=self.__time
294     if name == "set_up_time":
295     self.__set_up_time=float(value)
296     self.__cum_set_up_time+=self.__set_up_time
297 lgao 2596 if name == "net_time":
298     self.__net_time=float(value)
299     self.__cum_net_time+=self.__net_time
300 gross 2470 if name == "residual_norm":
301     self.__residual_norm=float(value)
302     if name == "converged":
303     self.__converged = (value == True)
304     def getDiagnostics(self, name):
305     """
306     Returns the diagnostic information C{name}
307    
308     @param name: name of diagnostic information where
309     - "num_iter": the number of iteration steps
310     - "cum_num_iter": the cumulative number of iteration steps
311     - "num_level": the number of level in multi level solver
312     - "num_inner_iter": the number of inner iteration steps
313     - "cum_num_inner_iter": the cumulative number of inner iteration steps
314     - "time": execution time
315     - "cum_time": cumulative execution time
316     - "set_up_time": time to set up of the solver, typically this includes factorization and reordering
317     - "cum_set_up_time": cumulative time to set up of the solver
318 lgao 2596 - "net_time": net execution time, excluding setup time for the solver and execution time for preconditioner
319     - "cum_net_time": cumulative net execution time
320 gross 2470 - "residual_norm": norm of the final residual
321     - "converged": return self.__converged
322 lgao 2596 @type name: C{str} in the list "num_iter", "num_level", "num_inner_iter", "time", "set_up_time", "net_time", "residual_norm", "converged".
323 gross 2470 @return: requested value. C{None} is returned if the value is undefined.
324     @note: If the solver has thrown an exception diagnostic values have an undefined status.
325     """
326     if name == "num_iter": return self.__num_iter
327     elif name == "cum_num_iter": return self.__cum_num_iter
328     elif name == "num_level": return self.__num_level
329     elif name == "num_inner_iter": return self.__num_inner_iter
330     elif name == "cum_num_inner_iter": return self.__cum_num_inner_iter
331     elif name == "time": return self.__time
332     elif name == "cum_time": return self.__cum_time
333     elif name == "set_up_time": return self.__set_up_time
334     elif name == "cum_set_up_time": return self.__cum_set_up_time
335 lgao 2596 elif name == "net_time": return self.__net_time
336     elif name == "cum_net_time": return self.__cum_net_time
337 gross 2470 elif name == "residual_norm": return self.__residual_norm
338     elif name == "converged": return self.__converged
339     else:
340     raise ValueError,"unknown diagnostic item %s"%name
341     def hasConverged(self):
342     """
343     Returns C{True} if the last solver call has been finalized successfully.
344     @note: if an exception has been thrown by the solver the status of this flag is undefined.
345     """
346     return self.getDiagnostics("converged")
347     def setCoarsening(self,method=0):
348     """
349     Sets the key of the coarsening method to be applied in AMG.
350    
351     @param method: selects the coarsening method .
352 gross 2474 @type method: in {SolverOptions.DEFAULT}, L{SolverOptions.YAIR_SHAPIRA_COARSENING},
353 gross 2470 L{SolverOptions.RUGE_STUEBEN_COARSENING}, L{SolverOptions.AGGREGATION_COARSENING}
354     """
355 gross 2474 if method==None: method=0
356 gross 2470 if not method in [self.DEFAULT, self.YAIR_SHAPIRA_COARSENING, self.RUGE_STUEBEN_COARSENING, self.AGGREGATION_COARSENING]:
357     raise ValueError,"unknown coarsening method %s"%method
358     self.__coarsening=method
359 artak 2524
360 gross 2470 def getCoarsening(self):
361     """
362     Returns the key of the coarsening algorithm to be applied AMG.
363    
364     @rtype: in the list L{SolverOptions.DEFAULT}, L{SolverOptions.YAIR_SHAPIRA_COARSENING},
365     L{SolverOptions.RUGE_STUEBEN_COARSENING}, L{SolverOptions.AGGREGATION_COARSENING}
366     """
367     return self.__coarsening
368 artak 2524
369     def setMinCoarseMatrixSize(self,size=500):
370     """
371     Sets the minumum size of the coarsest level matrix in AMG.
372    
373     @param size: minumum size of the coarsest level matrix .
374     @type size: positive C{int} or C{None}
375     """
376     size=int(size)
377     if size<0:
378     raise ValueError,"minumum size of the coarsest level matrix must be non-negative."
379     if size==None: size=500
380     self.__MinCoarseMatrixSize=size
381    
382     def getMinCoarseMatrixSize(self):
383     """
384     Returns the minumum size of the coarsest level matrix in AMG.
385    
386     @rtype: C{int}
387     """
388     return self.__MinCoarseMatrixSize
389    
390 gross 2470 def setPreconditioner(self, preconditioner=10):
391     """
392     Sets the preconditioner to be used.
393    
394     @param preconditioner: key of the preconditioner to be used.
395     @type preconditioner: in L{SolverOptions.SSOR}, L{SolverOptions.ILU0}, L{SolverOptions.ILUT}, L{SolverOptions.JACOBI},
396     L{SolverOptions.AMG}, L{SolverOptions.REC_ILU}, L{SolverOptions.GAUSS_SEIDEL}, L{SolverOptions.RILU},
397     L{SolverOptions.NO_PRECONDITIONER}
398     @note: Not all packages support all preconditioner. It can be assumed that a package makes a reasonable choice if it encounters
399     an unknown preconditioner.
400     """
401 gross 2474 if preconditioner==None: preconditioner=10
402 gross 2470 if not preconditioner in [ SolverOptions.SSOR, SolverOptions.ILU0, SolverOptions.ILUT, SolverOptions.JACOBI,
403     SolverOptions.AMG, SolverOptions.REC_ILU, SolverOptions.GAUSS_SEIDEL, SolverOptions.RILU,
404     SolverOptions.NO_PRECONDITIONER] :
405     raise ValueError,"unknown preconditioner %s"%preconditioner
406     self.__preconditioner=preconditioner
407     def getPreconditioner(self):
408     """
409     Returns key of the preconditioner to be used.
410    
411     @rtype: in the list L{SolverOptions.SSOR}, L{SolverOptions.ILU0}, L{SolverOptions.ILUT}, L{SolverOptions.JACOBI},
412     L{SolverOptions.AMG}, L{SolverOptions.REC_ILU}, L{SolverOptions.GAUSS_SEIDEL}, L{SolverOptions.RILU},
413     L{SolverOptions.NO_PRECONDITIONER}
414     """
415     return self.__preconditioner
416     def setSolverMethod(self, method=0):
417     """
418     Sets the solver method to be used. Use C{method}=C{SolverOptions.DIRECT} to indicate that a direct rather than an iterative
419     solver should be used and Use C{method}=C{SolverOptions.ITERATIVE} to indicate that an iterative rather than a direct
420     solver should be used.
421    
422     @param method: key of the solver method to be used.
423     @type method: in L{SolverOptions.DEFAULT}, L{SolverOptions.DIRECT}, L{SolverOptions.CHOLEVSKY}, L{SolverOptions.PCG},
424     L{SolverOptions.CR}, L{SolverOptions.CGS}, L{SolverOptions.BICGSTAB}, L{SolverOptions.SSOR},
425     L{SolverOptions.GMRES}, L{SolverOptions.PRES20}, L{SolverOptions.LUMPING}, L{SolverOptions.ITERATIVE},
426     L{SolverOptions.AMG}, L{SolverOptions.NONLINEAR_GMRES}, L{SolverOptions.TFQMR}, L{SolverOptions.MINRES},
427     L{SolverOptions.GAUSS_SEIDEL}
428     @note: Not all packages support all solvers. It can be assumed that a package makes a reasonable choice if it encounters
429     an unknown solver method.
430     """
431 gross 2474 if method==None: method=0
432 gross 2470 if not method in [ SolverOptions.DEFAULT, SolverOptions.DIRECT, SolverOptions.CHOLEVSKY, SolverOptions.PCG,
433     SolverOptions.CR, SolverOptions.CGS, SolverOptions.BICGSTAB, SolverOptions.SSOR,
434     SolverOptions.GMRES, SolverOptions.PRES20, SolverOptions.LUMPING, SolverOptions.ITERATIVE, SolverOptions.AMG,
435     SolverOptions.NONLINEAR_GMRES, SolverOptions.TFQMR, SolverOptions.MINRES, SolverOptions.GAUSS_SEIDEL]:
436     raise ValueError,"unknown solver method %s"%method
437     self.__method=method
438     def getSolverMethod(self):
439     """
440     Returns key of the solver method to be used.
441    
442     @rtype: in the list L{SolverOptions.DEFAULT}, L{SolverOptions.DIRECT}, L{SolverOptions.CHOLEVSKY}, L{SolverOptions.PCG},
443     L{SolverOptions.CR}, L{SolverOptions.CGS}, L{SolverOptions.BICGSTAB}, L{SolverOptions.SSOR},
444     L{SolverOptions.GMRES}, L{SolverOptions.PRES20}, L{SolverOptions.LUMPING}, L{SolverOptions.ITERATIVE},
445     L{SolverOptions.AMG}, L{SolverOptions.NONLINEAR_GMRES}, L{SolverOptions.TFQMR}, L{SolverOptions.MINRES},
446     L{SolverOptions.GAUSS_SEIDEL}
447     """
448     return self.__method
449    
450     def setPackage(self, package=0):
451     """
452     Sets the solver package to be used as a solver.
453    
454     @param package: key of the solver package to be used.
455     @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}
456     @note: Not all packages are support on all implementation. An exception may be thrown on some platforms if a particular is requested.
457     """
458 gross 2474 if package==None: package=0
459 gross 2470 if not package in [SolverOptions.DEFAULT, SolverOptions.PASO, SolverOptions.SUPER_LU, SolverOptions.PASTIX, SolverOptions.MKL, SolverOptions.UMFPACK, SolverOptions.TRILINOS]:
460     raise ValueError,"unknown solver package %s"%package
461     self.__package=package
462     def getPackage(self):
463     """
464     Returns the solver package key
465    
466     @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}
467     """
468     return self.__package
469     def setReordering(self,ordering=30):
470     """
471     Sets the key of the reordering method to be applied if supported by the solver. Some direct solvers support reordering
472     to optimize compute time and storage use during elimination.
473    
474     @param ordering: selects the reordering strategy.
475     @type ordering: in L{SolverOptions.NO_REORDERING}, L{SolverOptions.NO_REORDERING},
476     L{SolverOptions.NO_REORDERING}, L{SolverOptions.DEFAULT_REORDERING}
477     """
478     if not ordering in [self.NO_REORDERING, self.MINIMUM_FILL_IN, self.NESTED_DISSECTION, self.DEFAULT_REORDERING]:
479     raise ValueError,"unknown reordering strategy %s"%ordering
480     self.__reordering=ordering
481     def getReordering(self):
482     """
483     Returns the key of the reordering method to be applied if supported by the solver.
484    
485     @rtype: in the list L{SolverOptions.NO_REORDERING}, L{SolverOptions.NO_REORDERING},
486     L{SolverOptions.NO_REORDERING}, L{SolverOptions.DEFAULT_REORDERING}
487     """
488     return self.__reordering
489     def setRestart(self,restart=None):
490     """
491     Sets the number of iterations steps after which GMRES is performing a restart.
492    
493     @param restart: number of iteration steps after which to perform a restart. If equal to C{None} no
494     restart is performed.
495     @type restart: C{int} or C{None}
496     """
497     if restart == None:
498     self.__restart=restart
499     else:
500     restart=int(restart)
501     if restart<1:
502     raise ValueError,"restart must be positive."
503     self.__restart=restart
504 gross 2474
505 gross 2470 def getRestart(self):
506     """
507     Returns the number of iterations steps after which GMRES is performing a restart.
508     If C{None} is returned no restart is performed.
509    
510     @rtype: C{int} or C{None}
511     """
512     if self.__restart < 0:
513     return None
514     else:
515     return self.__restart
516 gross 2474 def _getRestartForC(self):
517     r=self.getRestart()
518     if r == None:
519     return -1
520     else:
521     return r
522 gross 2470 def setTruncation(self,truncation=20):
523     """
524     Sets the number of residuals in GMRES to be stored for orthogonalization. The more residuals are stored
525     the faster GMRES converged but
526    
527     @param truncation: truncation
528     @type truncation: C{int}
529     """
530     truncation=int(truncation)
531     if truncation<1:
532     raise ValueError,"truncation must be positive."
533     self.__truncation=truncation
534     def getTruncation(self):
535     """
536     Returns the number of residuals in GMRES to be stored for orthogonalization
537    
538     @rtype: C{int}
539     """
540     return self.__truncation
541     def setInnerIterMax(self,iter_max=10):
542     """
543     Sets the maximum number of iteration steps for the inner iteration.
544    
545     @param iter_max: maximum number of inner iterations
546     @type iter_max: C{int}
547     """
548     iter_max=int(iter_max)
549     if iter_max<1:
550     raise ValueError,"maximum number of inner iteration must be positive."
551     self.__inner_iter_max=iter_max
552     def getInnerIterMax(self):
553     """
554     Returns maximum number of inner iteration steps
555    
556     @rtype: C{int}
557     """
558     return self.__inner_iter_max
559 gross 2479 def setIterMax(self,iter_max=100000):
560 gross 2470 """
561     Sets the maximum number of iteration steps
562    
563     @param iter_max: maximum number of iteration steps
564     @type iter_max: C{int}
565     """
566     iter_max=int(iter_max)
567     if iter_max<1:
568     raise ValueError,"maximum number of iteration steps must be positive."
569     self.__iter_max=iter_max
570     def getIterMax(self):
571     """
572     Returns maximum number of iteration steps
573    
574     @rtype: C{int}
575     """
576     return self.__iter_max
577     def setLevelMax(self,level_max=10):
578     """
579     Sets the maximum number of coarsening levels to be used in an algebraic multi level solver or preconditioner
580    
581     @param level_max: maximum number of levels
582     @type level_max: C{int}
583     """
584     level_max=int(level_max)
585     if level_max<0:
586     raise ValueError,"maximum number of coarsening levels must be non-negative."
587     self.__level_max=level_max
588     def getLevelMax(self):
589     """
590     Returns the maximum number of coarsening levels to be used in an algebraic multi level solver or preconditioner
591    
592     @rtype: C{int}
593     """
594     return self.__level_max
595     def setCoarseningThreshold(self,theta=0.05):
596     """
597     Sets the threshold for coarsening in the algebraic multi level solver or preconditioner
598    
599     @param theta: threshold for coarsening
600     @type theta: positive C{float}
601     """
602     theta=float(theta)
603     if theta<0 or theta>1:
604     raise ValueError,"threshold must be non-negative and less or equal 1."
605     self.__coarsening_threshold=theta
606     def getCoarseningThreshold(self):
607     """
608     Returns the threshold for coarsening in the algebraic multi level solver or preconditioner
609    
610     @rtype: C{float}
611     """
612     return self.__coarsening_threshold
613     def setNumSweeps(self,sweeps=2):
614     """
615     Sets the number of sweeps in a Jacobi or Gauss-Seidel/SOR preconditioner.
616    
617     @param sweeps: number of sweeps
618     @type theta: positive C{int}
619     """
620     sweeps=int(sweeps)
621     if sweeps<1:
622     raise ValueError,"number of sweeps must be positive."
623     self.__sweeps=sweeps
624     def getNumSweeps(self):
625     """
626     Returns the number of sweeps in a Jacobi or Gauss-Seidel/SOR preconditioner.
627    
628     @rtype: C{int}
629     """
630     return self.__sweeps
631     def setNumPreSweeps(self,sweeps=2):
632     """
633     Sets the number of sweeps in the pre-smoothing step of a multi level solver or preconditioner
634    
635     @param sweeps: number of sweeps
636     @type theta: positive C{int}
637     """
638     sweeps=int(sweeps)
639     if sweeps<1:
640     raise ValueError,"number of sweeps must be positive."
641     self.__pre_sweeps=sweeps
642     def getNumPreSweeps(self):
643     """
644     Returns he number of sweeps in the pre-smoothing step of a multi level solver or preconditioner
645    
646     @rtype: C{int}
647     """
648     return self.__pre_sweeps
649     def setNumPostSweeps(self,sweeps=2):
650     """
651     Sets the number of sweeps in the post-smoothing step of a multi level solver or preconditioner
652    
653     @param sweeps: number of sweeps
654     @type theta: positive C{int}
655     """
656     sweeps=int(sweeps)
657     if sweeps<1:
658     raise ValueError,"number of sweeps must be positive."
659     self.__post_sweeps=sweeps
660     def getNumPostSweeps(self):
661     """
662     Returns he number of sweeps in the post-smoothing step of a multi level solver or preconditioner
663    
664     @rtype: C{int}
665     """
666     return self.__post_sweeps
667    
668     def setTolerance(self,rtol=1.e-8):
669     """
670     Sets the relative tolerance for the solver
671    
672     @param rtol: relative tolerance
673     @type rtol: non-negative C{float}
674     """
675     rtol=float(rtol)
676     if rtol<0 or rtol>1:
677     raise ValueError,"tolerance must be non-negative and less or equal 1."
678     self.__tolerance=rtol
679     def getTolerance(self):
680     """
681     Returns the relative tolerance for the solver
682    
683     @rtype: C{float}
684     """
685     return self.__tolerance
686     def setAbsoluteTolerance(self,atol=0.):
687     """
688     Sets the absolute tolerance for the solver
689    
690     @param atol: absolute tolerance
691     @type atol: non-negative C{float}
692     """
693     atol=float(atol)
694     if atol<0:
695     raise ValueError,"tolerance must be non-negative."
696     self.__absolute_tolerance=atol
697     def getAbsoluteTolerance(self):
698     """
699     Returns the absolute tolerance for the solver
700    
701     @rtype: C{float}
702     """
703     return self.__absolute_tolerance
704    
705     def setInnerTolerance(self,rtol=0.9):
706     """
707     Sets the relative tolerance for an inner iteration scheme for instance
708     on the coarsest level in a multi-level scheme.
709    
710     @param rtol: inner relative tolerance
711     @type rtol: positive C{float}
712     """
713     rtol=float(rtol)
714     if rtol<=0 or rtol>1:
715     raise ValueError,"tolerance must be positive and less or equal 1."
716     self.__inner_tolerance=rtol
717     def getInnerTolerance(self):
718     """
719     Returns the relative tolerance for an inner iteration scheme
720    
721     @rtype: C{float}
722     """
723     return self.__inner_tolerance
724     def setDropTolerance(self,drop_tol=0.01):
725     """
726     Sets the relative drop tolerance in ILUT
727    
728     @param drop_tol: drop tolerance
729     @type drop_tol: positive C{float}
730     """
731     drop_tol=float(drop_tol)
732     if drop_tol<=0 or drop_tol>1:
733     raise ValueError,"drop tolerance must be positive and less or equal 1."
734     self.__drop_tolerance=drop_tol
735     def getDropTolerance(self):
736     """
737     Returns the relative drop tolerance in ILUT
738    
739     @rtype: C{float}
740     """
741     return self.__drop_tolerance
742     def setDropStorage(self,storage=2.):
743     """
744     Sets the maximum allowed increase in storage for ILUT. C{storage}=2 would mean that
745     a doubling of the storage needed for the coefficient matrix is allowed in the ILUT factorization.
746    
747     @param storage: allowed storage increase
748     @type storage: C{float}
749     """
750     storage=float(storage)
751     if storage<1:
752     raise ValueError,"allowed storage increase must be greater or equal to 1."
753     self.__drop_storage=storage
754     def getDropStorage(self):
755    
756     """
757     Returns the maximum allowed increase in storage for ILUT
758    
759     @rtype: C{float}
760     """
761     return self.__drop_storage
762     def setRelaxationFactor(self,factor=0.3):
763     """
764     Sets the relaxation factor used to add dropped elements in RILU to the main diagonal.
765    
766     @param factor: relaxation factor
767     @type factor: C{float}
768     @note: RILU with a relaxation factor 0 is identical to ILU0
769     """
770     factor=float(factor)
771     if factor<0:
772     raise ValueError,"relaxation factor must be non-negative."
773     self.__relaxation=factor
774     def getRelaxationFactor(self):
775    
776     """
777     Returns the relaxation factor used to add dropped elements in RILU to the main diagonal.
778    
779     @rtype: C{float}
780     """
781     return self.__relaxation
782     def isSymmetric(self):
783     """
784     Checks if symmetry of the coefficient matrix is indicated.
785    
786     @return: True if a symmetric PDE is indicated, False otherwise
787     @rtype: C{bool}
788     """
789     return self.__symmetric
790     def setSymmetryOn(self):
791     """
792     Sets the symmetry flag to indicate that the coefficient matrix is symmetric.
793     """
794     self.__symmetric=True
795     def setSymmetryOff(self):
796     """
797     Clears the symmetry flag for the coefficient matrix.
798     """
799     self.__symmetric=False
800     def setSymmetry(self,flag=False):
801     """
802     Sets the symmetry flag for the coefficient matrix to C{flag}.
803    
804     @param flag: If True, the symmetry flag is set otherwise reset.
805     @type flag: C{bool}
806     """
807     if flag:
808     self.setSymmetryOn()
809     else:
810     self.setSymmetryOff()
811     def isVerbose(self):
812     """
813     Returns C{True} if the solver is expected to be verbose.
814    
815     @return: True if verbosity of switched on.
816     @rtype: C{bool}
817     """
818     return self.__verbose
819    
820     def setVerbosityOn(self):
821     """
822     Switches the verbosity of the solver on.
823     """
824     self.__verbose=True
825     def setVerbosityOff(self):
826     """
827     Switches the verbosity of the solver off.
828     """
829     self.__verbose=False
830 gross 2474 def setVerbosity(self,verbose=False):
831 gross 2470 """
832     Sets the verbosity flag for the solver to C{flag}.
833    
834     @param flag: If C{True}, the verbosity of the solver is switched on.
835     @type flag: C{bool}
836     """
837     if verbose:
838     self.setVerbosityOn()
839     else:
840     self.setVerbosityOff()
841 gross 2474
842 gross 2470 def adaptInnerTolerance(self):
843     """
844     Returns C{True} if the tolerance of the inner solver is selected automatically.
845     Otherwise the inner tolerance set by L{setInnerTolerance} is used.
846    
847     @return: C{True} if inner tolerance adaption is chosen.
848     @rtype: C{bool}
849     """
850     return self.__adapt_inner_tolerance
851    
852     def setInnerToleranceAdaptionOn(self):
853     """
854     Switches the automatic selection of inner tolerance on
855     """
856     self.__adapt_inner_tolerance=True
857     def setInnerToleranceAdaptionOff(self):
858     """
859     Switches the automatic selection of inner tolerance off.
860     """
861     self.__adapt_inner_tolerance=False
862     def setInnerToleranceAdaption(self,adapt=True):
863     """
864     Sets a flag to indicate automatic selection of the inner tolerance.
865    
866     @param adapt: If C{True}, the inner tolerance is selected automatically.
867     @type adapt: C{bool}
868     """
869     if adapt:
870     self.setInnerToleranceAdaptionOn()
871     else:
872     self.setInnerToleranceAdaptionOff()
873    
874     def acceptConvergenceFailure(self):
875     """
876     Returns C{True} if a failure to meet the stopping criteria within the
877     given number of iteration steps is not raising in exception. This is useful
878     if a solver is used in a non-linear context where the non-linear solver can
879     continue even if the returned the solution does not necessarily meet the
880     stopping criteria. One can use the L{hasConverged} method to check if the
881     last call to the solver was successful.
882    
883     @return: C{True} if a failure to achieve convergence is accepted.
884     @rtype: C{bool}
885     """
886     return self.__accept_convergence_failure
887    
888     def setAcceptanceConvergenceFailureOn(self):
889     """
890     Switches the acceptance of a failure of convergence on
891     """
892     self.__accept_convergence_failure=True
893     def setAcceptanceConvergenceFailureOff(self):
894     """
895     Switches the acceptance of a failure of convergence off.
896     """
897     self.__accept_convergence_failure=False
898     def setAcceptanceConvergenceFailure(self,accept=False):
899     """
900     Sets a flag to indicate the acceptance of a failure of convergence.
901    
902     @param accept: If C{True}, any failure to achieve convergence is accepted.
903     @type accept: C{bool}
904     """
905     if accept:
906     self.setAcceptanceConvergenceFailureOn()
907     else:
908     self.setAcceptanceConvergenceFailureOff()
909    
910 jgs 148 class IllegalCoefficient(ValueError):
911 jgs 102 """
912 caltinay 2169 Exception that is raised if an illegal coefficient of the general or
913     particular PDE is requested.
914 jgs 148 """
915 gross 1072 pass
916 jgs 102
917 jgs 148 class IllegalCoefficientValue(ValueError):
918 jgs 102 """
919 caltinay 2169 Exception that is raised if an incorrect value for a coefficient is used.
920 jgs 148 """
921 gross 1072 pass
922 jgs 122
923 gross 1072 class IllegalCoefficientFunctionSpace(ValueError):
924     """
925 caltinay 2169 Exception that is raised if an incorrect function space for a coefficient
926     is used.
927 gross 1072 """
928    
929 jgs 148 class UndefinedPDEError(ValueError):
930     """
931 caltinay 2169 Exception that is raised if a PDE is not fully defined yet.
932 jgs 148 """
933 gross 1072 pass
934 jgs 102
935 gross 1841 class PDECoef(object):
936 jgs 102 """
937 caltinay 2169 A class for describing a PDE coefficient.
938 jgs 149
939 caltinay 2169 @cvar INTERIOR: indicator that coefficient is defined on the interior of
940     the PDE domain
941     @cvar BOUNDARY: indicator that coefficient is defined on the boundary of
942     the PDE domain
943     @cvar CONTACT: indicator that coefficient is defined on the contact region
944     within the PDE domain
945     @cvar INTERIOR_REDUCED: indicator that coefficient is defined on the
946     interior of the PDE domain using a reduced
947     integration order
948     @cvar BOUNDARY_REDUCED: indicator that coefficient is defined on the
949     boundary of the PDE domain using a reduced
950     integration order
951     @cvar CONTACT_REDUCED: indicator that coefficient is defined on the contact
952     region within the PDE domain using a reduced
953     integration order
954     @cvar SOLUTION: indicator that coefficient is defined through a solution of
955     the PDE
956     @cvar REDUCED: indicator that coefficient is defined through a reduced
957     solution of the PDE
958     @cvar BY_EQUATION: indicator that the dimension of the coefficient shape is
959     defined by the number of PDE equations
960     @cvar BY_SOLUTION: indicator that the dimension of the coefficient shape is
961     defined by the number of PDE solutions
962     @cvar BY_DIM: indicator that the dimension of the coefficient shape is
963     defined by the spatial dimension
964     @cvar OPERATOR: indicator that the the coefficient alters the operator of
965     the PDE
966     @cvar RIGHTHANDSIDE: indicator that the the coefficient alters the right
967     hand side of the PDE
968     @cvar BOTH: indicator that the the coefficient alters the operator as well
969     as the right hand side of the PDE
970 jgs 149
971 jgs 102 """
972     INTERIOR=0
973     BOUNDARY=1
974     CONTACT=2
975 jgs 149 SOLUTION=3
976 jgs 150 REDUCED=4
977 jgs 149 BY_EQUATION=5
978     BY_SOLUTION=6
979     BY_DIM=7
980     OPERATOR=10
981     RIGHTHANDSIDE=11
982     BOTH=12
983 gross 1072 INTERIOR_REDUCED=13
984     BOUNDARY_REDUCED=14
985     CONTACT_REDUCED=15
986 jgs 149
987 gross 1072 def __init__(self, where, pattern, altering):
988 jgs 102 """
989 caltinay 2169 Initialises a PDE coefficient type.
990 jgs 151
991 jgs 149 @param where: describes where the coefficient lives
992 caltinay 2169 @type where: one of L{INTERIOR}, L{BOUNDARY}, L{CONTACT}, L{SOLUTION},
993     L{REDUCED}, L{INTERIOR_REDUCED}, L{BOUNDARY_REDUCED},
994     L{CONTACT_REDUCED}
995     @param pattern: describes the shape of the coefficient and how the shape
996     is built for a given spatial dimension and numbers of
997     equations and solutions in then PDE. For instance,
998     (L{BY_EQUATION},L{BY_SOLUTION},L{BY_DIM}) describes a
999     rank 3 coefficient which is instantiated as shape (3,2,2)
1000     in case of three equations and two solution components
1001     on a 2-dimensional domain. In the case of single equation
1002     and a single solution component the shape components
1003     marked by L{BY_EQUATION} or L{BY_SOLUTION} are dropped.
1004     In this case the example would be read as (2,).
1005 jgs 149 @type pattern: C{tuple} of L{BY_EQUATION}, L{BY_SOLUTION}, L{BY_DIM}
1006 caltinay 2169 @param altering: indicates what part of the PDE is altered if the
1007     coefficient is altered
1008 jgs 149 @type altering: one of L{OPERATOR}, L{RIGHTHANDSIDE}, L{BOTH}
1009 jgs 102 """
1010 gross 1841 super(PDECoef, self).__init__()
1011 jgs 102 self.what=where
1012     self.pattern=pattern
1013     self.altering=altering
1014 jgs 108 self.resetValue()
1015 jgs 102
1016 jgs 108 def resetValue(self):
1017     """
1018 caltinay 2169 Resets the coefficient value to the default.
1019 jgs 108 """
1020     self.value=escript.Data()
1021    
1022 jgs 149 def getFunctionSpace(self,domain,reducedEquationOrder=False,reducedSolutionOrder=False):
1023 jgs 102 """
1024 caltinay 2169 Returns the L{FunctionSpace<escript.FunctionSpace>} of the coefficient.
1025 jgs 102
1026 jgs 149 @param domain: domain on which the PDE uses the coefficient
1027     @type domain: L{Domain<escript.Domain>}
1028 caltinay 2169 @param reducedEquationOrder: True to indicate that reduced order is used
1029     to represent the equation
1030 gross 720 @type reducedEquationOrder: C{bool}
1031 caltinay 2169 @param reducedSolutionOrder: True to indicate that reduced order is used
1032     to represent the solution
1033 gross 720 @type reducedSolutionOrder: C{bool}
1034 caltinay 2169 @return: L{FunctionSpace<escript.FunctionSpace>} of the coefficient
1035     @rtype: L{FunctionSpace<escript.FunctionSpace>}
1036 jgs 102 """
1037 jgs 151 if self.what==self.INTERIOR:
1038 jgs 149 return escript.Function(domain)
1039 gross 1072 elif self.what==self.INTERIOR_REDUCED:
1040     return escript.ReducedFunction(domain)
1041 jgs 151 elif self.what==self.BOUNDARY:
1042 jgs 149 return escript.FunctionOnBoundary(domain)
1043 gross 1072 elif self.what==self.BOUNDARY_REDUCED:
1044     return escript.ReducedFunctionOnBoundary(domain)
1045 jgs 151 elif self.what==self.CONTACT:
1046 jgs 149 return escript.FunctionOnContactZero(domain)
1047 gross 1072 elif self.what==self.CONTACT_REDUCED:
1048     return escript.ReducedFunctionOnContactZero(domain)
1049 jgs 151 elif self.what==self.SOLUTION:
1050 jgs 149 if reducedEquationOrder and reducedSolutionOrder:
1051     return escript.ReducedSolution(domain)
1052     else:
1053     return escript.Solution(domain)
1054 jgs 151 elif self.what==self.REDUCED:
1055     return escript.ReducedSolution(domain)
1056 jgs 102
1057 jgs 108 def getValue(self):
1058     """
1059 caltinay 2169 Returns the value of the coefficient.
1060 jgs 149
1061 caltinay 2169 @return: value of the coefficient
1062     @rtype: L{Data<escript.Data>}
1063 jgs 108 """
1064     return self.value
1065 jgs 148
1066 jgs 149 def setValue(self,domain,numEquations=1,numSolutions=1,reducedEquationOrder=False,reducedSolutionOrder=False,newValue=None):
1067 jgs 108 """
1068 caltinay 2169 Sets the value of the coefficient to a new value.
1069 jgs 149
1070     @param domain: domain on which the PDE uses the coefficient
1071     @type domain: L{Domain<escript.Domain>}
1072     @param numEquations: number of equations of the PDE
1073     @type numEquations: C{int}
1074     @param numSolutions: number of components of the PDE solution
1075     @type numSolutions: C{int}
1076 caltinay 2169 @param reducedEquationOrder: True to indicate that reduced order is used
1077     to represent the equation
1078 gross 720 @type reducedEquationOrder: C{bool}
1079 caltinay 2169 @param reducedSolutionOrder: True to indicate that reduced order is used
1080     to represent the solution
1081 gross 720 @type reducedSolutionOrder: C{bool}
1082 jgs 149 @param newValue: number of components of the PDE solution
1083 caltinay 2169 @type newValue: any object that can be converted into a
1084     L{Data<escript.Data>} object with the appropriate shape
1085     and L{FunctionSpace<escript.FunctionSpace>}
1086     @raise IllegalCoefficientValue: if the shape of the assigned value does
1087     not match the shape of the coefficient
1088     @raise IllegalCoefficientFunctionSpace: if unable to interpolate value
1089     to appropriate function space
1090 jgs 108 """
1091 jgs 148 if newValue==None:
1092     newValue=escript.Data()
1093     elif isinstance(newValue,escript.Data):
1094     if not newValue.isEmpty():
1095 gross 1072 if not newValue.getFunctionSpace() == self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder):
1096     try:
1097     newValue=escript.Data(newValue,self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder))
1098     except:
1099     raise IllegalCoefficientFunctionSpace,"Unable to interpolate coefficient to function space %s"%self.getFunctionSpace(domain)
1100 jgs 148 else:
1101 jgs 149 newValue=escript.Data(newValue,self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder))
1102 jgs 148 if not newValue.isEmpty():
1103     if not self.getShape(domain,numEquations,numSolutions)==newValue.getShape():
1104 jgs 149 raise IllegalCoefficientValue,"Expected shape of coefficient is %s but actual shape is %s."%(self.getShape(domain,numEquations,numSolutions),newValue.getShape())
1105 jgs 108 self.value=newValue
1106 jgs 148
1107 jgs 102 def isAlteringOperator(self):
1108     """
1109 caltinay 2169 Checks if the coefficient alters the operator of the PDE.
1110 jgs 149
1111 caltinay 2169 @return: True if the operator of the PDE is changed when the
1112     coefficient is changed
1113     @rtype: C{bool}
1114     """
1115 jgs 102 if self.altering==self.OPERATOR or self.altering==self.BOTH:
1116     return not None
1117     else:
1118     return None
1119    
1120     def isAlteringRightHandSide(self):
1121     """
1122 caltinay 2169 Checks if the coefficient alters the right hand side of the PDE.
1123 jgs 149
1124 caltinay 2169 @rtype: C{bool}
1125     @return: True if the right hand side of the PDE is changed when the
1126     coefficient is changed, C{None} otherwise.
1127     """
1128 jgs 102 if self.altering==self.RIGHTHANDSIDE or self.altering==self.BOTH:
1129     return not None
1130     else:
1131     return None
1132    
1133 jgs 148 def estimateNumEquationsAndNumSolutions(self,domain,shape=()):
1134 jgs 102 """
1135 caltinay 2169 Tries to estimate the number of equations and number of solutions if
1136     the coefficient has the given shape.
1137 jgs 102
1138 jgs 149 @param domain: domain on which the PDE uses the coefficient
1139     @type domain: L{Domain<escript.Domain>}
1140     @param shape: suggested shape of the coefficient
1141     @type shape: C{tuple} of C{int} values
1142 caltinay 2169 @return: the number of equations and number of solutions of the PDE if
1143     the coefficient has given shape. If no appropriate numbers
1144     could be identified, C{None} is returned
1145 jgs 149 @rtype: C{tuple} of two C{int} values or C{None}
1146 jgs 102 """
1147 jgs 148 dim=domain.getDim()
1148 jgs 102 if len(shape)>0:
1149     num=max(shape)+1
1150     else:
1151     num=1
1152     search=[]
1153 jgs 149 if self.definesNumEquation() and self.definesNumSolutions():
1154     for u in range(num):
1155     for e in range(num):
1156     search.append((e,u))
1157     search.sort(self.__CompTuple2)
1158     for item in search:
1159 jgs 148 s=self.getShape(domain,item[0],item[1])
1160 jgs 102 if len(s)==0 and len(shape)==0:
1161     return (1,1)
1162     else:
1163     if s==shape: return item
1164 jgs 149 elif self.definesNumEquation():
1165     for e in range(num,0,-1):
1166     s=self.getShape(domain,e,0)
1167     if len(s)==0 and len(shape)==0:
1168     return (1,None)
1169     else:
1170     if s==shape: return (e,None)
1171    
1172     elif self.definesNumSolutions():
1173     for u in range(num,0,-1):
1174     s=self.getShape(domain,0,u)
1175     if len(s)==0 and len(shape)==0:
1176     return (None,1)
1177     else:
1178     if s==shape: return (None,u)
1179 jgs 102 return None
1180 caltinay 2169
1181 jgs 149 def definesNumSolutions(self):
1182     """
1183 caltinay 2169 Checks if the coefficient allows to estimate the number of solution
1184     components.
1185 jgs 102
1186 caltinay 2169 @return: True if the coefficient allows an estimate of the number of
1187     solution components, False otherwise
1188 jgs 149 @rtype: C{bool}
1189     """
1190     for i in self.pattern:
1191     if i==self.BY_SOLUTION: return True
1192     return False
1193    
1194     def definesNumEquation(self):
1195     """
1196 caltinay 2169 Checks if the coefficient allows to estimate the number of equations.
1197 jgs 149
1198 caltinay 2169 @return: True if the coefficient allows an estimate of the number of
1199     equations, False otherwise
1200 jgs 149 @rtype: C{bool}
1201     """
1202     for i in self.pattern:
1203     if i==self.BY_EQUATION: return True
1204     return False
1205    
1206     def __CompTuple2(self,t1,t2):
1207     """
1208 caltinay 2169 Compares two tuples of possible number of equations and number of
1209     solutions.
1210 jgs 149
1211 caltinay 2169 @param t1: the first tuple
1212     @param t2: the second tuple
1213     @return: 0, 1, or -1
1214 jgs 149 """
1215    
1216     dif=t1[0]+t1[1]-(t2[0]+t2[1])
1217     if dif<0: return 1
1218     elif dif>0: return -1
1219     else: return 0
1220    
1221 jgs 148 def getShape(self,domain,numEquations=1,numSolutions=1):
1222 jgs 149 """
1223 caltinay 2169 Builds the required shape of the coefficient.
1224 jgs 102
1225 jgs 149 @param domain: domain on which the PDE uses the coefficient
1226     @type domain: L{Domain<escript.Domain>}
1227     @param numEquations: number of equations of the PDE
1228     @type numEquations: C{int}
1229     @param numSolutions: number of components of the PDE solution
1230     @type numSolutions: C{int}
1231     @return: shape of the coefficient
1232     @rtype: C{tuple} of C{int} values
1233     """
1234     dim=domain.getDim()
1235     s=()
1236     for i in self.pattern:
1237     if i==self.BY_EQUATION:
1238 jgs 148 if numEquations>1: s=s+(numEquations,)
1239 jgs 149 elif i==self.BY_SOLUTION:
1240 jgs 148 if numSolutions>1: s=s+(numSolutions,)
1241 jgs 102 else:
1242     s=s+(dim,)
1243 jgs 149 return s
1244 jgs 102
1245 gross 2470 #====================================================================================================================
1246    
1247 gross 1841 class LinearProblem(object):
1248 jgs 102 """
1249 caltinay 2169 This is the base class to define a general linear PDE-type problem for
1250     for an unknown function M{u} on a given domain defined through a
1251     L{Domain<escript.Domain>} object. The problem can be given as a single
1252     equation or as a system of equations.
1253 jgs 102
1254 caltinay 2169 The class assumes that some sort of assembling process is required to form
1255     a problem of the form
1256 jgs 151
1257 caltinay 2169 M{L u=f}
1258 jgs 102
1259 caltinay 2169 where M{L} is an operator and M{f} is the right hand side. This operator
1260     problem will be solved to get the unknown M{u}.
1261    
1262 jgs 102 """
1263 jgs 148 def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):
1264 jgs 102 """
1265 caltinay 2169 Initializes a linear problem.
1266 jgs 102
1267 jgs 148 @param domain: domain of the PDE
1268 jgs 149 @type domain: L{Domain<escript.Domain>}
1269 caltinay 2169 @param numEquations: number of equations. If C{None} the number of
1270     equations is extracted from the coefficients.
1271     @param numSolutions: number of solution components. If C{None} the number
1272     of solution components is extracted from the
1273     coefficients.
1274     @param debug: if True debug information is printed
1275 jgs 148
1276 jgs 102 """
1277 gross 1841 super(LinearProblem, self).__init__()
1278 jgs 102
1279 jgs 148 self.__debug=debug
1280 jgs 104 self.__domain=domain
1281     self.__numEquations=numEquations
1282     self.__numSolutions=numSolutions
1283 gross 1841 self.__altered_coefficients=False
1284 jgs 149 self.__reduce_equation_order=False
1285     self.__reduce_solution_order=False
1286 gross 2474 self.__sym=False
1287     self.setSolverOptions()
1288 gross 1841 self.__is_RHS_valid=False
1289     self.__is_operator_valid=False
1290     self.__COEFFICIENTS={}
1291 gross 2474 self.__solution_rtol=1.e99
1292     self.__solution_atol=1.e99
1293     self.setSymmetryOff()
1294 gross 1841 # initialize things:
1295     self.resetAllCoefficients()
1296 gross 2474 self.initializeSystem()
1297 caltinay 2169 # ==========================================================================
1298 jgs 148 # general stuff:
1299 caltinay 2169 # ==========================================================================
1300 jgs 148 def __str__(self):
1301 jgs 149 """
1302 caltinay 2169 Returns a string representation of the PDE.
1303 jgs 149
1304     @return: a simple representation of the PDE
1305     @rtype: C{str}
1306     """
1307 gross 1841 return "<LinearProblem %d>"%id(self)
1308 caltinay 2169 # ==========================================================================
1309 jgs 148 # debug :
1310 caltinay 2169 # ==========================================================================
1311 jgs 148 def setDebugOn(self):
1312 jgs 149 """
1313 caltinay 2169 Switches debug output on.
1314 jgs 108 """
1315 jgs 148 self.__debug=not None
1316    
1317     def setDebugOff(self):
1318 jgs 108 """
1319 caltinay 2169 Switches debug output off.
1320 jgs 148 """
1321     self.__debug=None
1322 jgs 108
1323 gross 2100 def setDebug(self, flag):
1324     """
1325 caltinay 2169 Switches debug output on if C{flag} is True otherwise it is switched off.
1326    
1327 gross 2100 @param flag: desired debug status
1328     @type flag: C{bool}
1329     """
1330     if flag:
1331     self.setDebugOn()
1332     else:
1333     self.setDebugOff()
1334    
1335 jgs 148 def trace(self,text):
1336     """
1337 caltinay 2169 Prints the text message if debug mode is switched on.
1338    
1339     @param text: message to be printed
1340 jgs 149 @type text: C{string}
1341 jgs 102 """
1342 jgs 148 if self.__debug: print "%s: %s"%(str(self),text)
1343 jgs 102
1344 caltinay 2169 # ==========================================================================
1345 jgs 148 # some service functions:
1346 caltinay 2169 # ==========================================================================
1347 gross 1841 def introduceCoefficients(self,**coeff):
1348     """
1349 caltinay 2169 Introduces new coefficients into the problem.
1350 gross 1841
1351 caltinay 2169 Use::
1352 gross 1841
1353 caltinay 2169 p.introduceCoefficients(A=PDECoef(...), B=PDECoef(...))
1354 gross 1841
1355     to introduce the coefficients M{A} ans M{B}.
1356     """
1357     for name, type in coeff.items():
1358     if not isinstance(type,PDECoef):
1359     raise ValueError,"coefficient %s has no type."%name
1360     self.__COEFFICIENTS[name]=type
1361     self.__COEFFICIENTS[name].resetValue()
1362     self.trace("coefficient %s has been introduced."%name)
1363    
1364 jgs 148 def getDomain(self):
1365 jgs 102 """
1366 caltinay 2169 Returns the domain of the PDE.
1367 jgs 102
1368 jgs 149 @return: the domain of the PDE
1369     @rtype: L{Domain<escript.Domain>}
1370 jgs 108 """
1371 jgs 148 return self.__domain
1372 gross 2535 def getDomainStatus(self):
1373     """
1374     Return the status indicator of the domain
1375     """
1376     return self.getDomain().getStatus()
1377 jgs 122
1378 gross 2535 def getSystemStatus(self):
1379     """
1380     Return the domain status used to build the current system
1381     """
1382     return self.__system_status
1383     def setSystemStatus(self,status=None):
1384     """
1385     Sets the system status to C{status} if C{status} is not present the
1386     current status of the domain is used.
1387     """
1388     if status == None:
1389     self.__system_status=self.getDomainStatus()
1390     else:
1391     self.__system_status=status
1392    
1393 jgs 148 def getDim(self):
1394 jgs 108 """
1395 caltinay 2169 Returns the spatial dimension of the PDE.
1396 jgs 108
1397 jgs 149 @return: the spatial dimension of the PDE domain
1398     @rtype: C{int}
1399 jgs 148 """
1400     return self.getDomain().getDim()
1401 jgs 102
1402 jgs 148 def getNumEquations(self):
1403     """
1404 caltinay 2169 Returns the number of equations.
1405 jgs 102
1406 jgs 149 @return: the number of equations
1407     @rtype: C{int}
1408 caltinay 2169 @raise UndefinedPDEError: if the number of equations is not specified yet
1409 jgs 148 """
1410     if self.__numEquations==None:
1411 gross 1859 if self.__numSolutions==None:
1412     raise UndefinedPDEError,"Number of equations is undefined. Please specify argument numEquations."
1413     else:
1414     self.__numEquations=self.__numSolutions
1415     return self.__numEquations
1416 jgs 147
1417 jgs 148 def getNumSolutions(self):
1418     """
1419 caltinay 2169 Returns the number of unknowns.
1420 jgs 147
1421 jgs 149 @return: the number of unknowns
1422     @rtype: C{int}
1423 caltinay 2169 @raise UndefinedPDEError: if the number of unknowns is not specified yet
1424 jgs 148 """
1425     if self.__numSolutions==None:
1426 gross 1859 if self.__numEquations==None:
1427     raise UndefinedPDEError,"Number of solution is undefined. Please specify argument numSolutions."
1428     else:
1429     self.__numSolutions=self.__numEquations
1430     return self.__numSolutions
1431 jgs 148
1432 jgs 149 def reduceEquationOrder(self):
1433     """
1434 caltinay 2169 Returns the status of order reduction for the equation.
1435 jgs 149
1436 caltinay 2169 @return: True if reduced interpolation order is used for the
1437     representation of the equation, False otherwise
1438 jgs 149 @rtype: L{bool}
1439     """
1440     return self.__reduce_equation_order
1441    
1442     def reduceSolutionOrder(self):
1443     """
1444 caltinay 2169 Returns the status of order reduction for the solution.
1445 jgs 149
1446 caltinay 2169 @return: True if reduced interpolation order is used for the
1447     representation of the solution, False otherwise
1448 jgs 149 @rtype: L{bool}
1449     """
1450     return self.__reduce_solution_order
1451 jgs 151
1452 jgs 108 def getFunctionSpaceForEquation(self):
1453     """
1454 caltinay 2169 Returns the L{FunctionSpace<escript.FunctionSpace>} used to discretize
1455     the equation.
1456 jgs 148
1457 jgs 149 @return: representation space of equation
1458     @rtype: L{FunctionSpace<escript.FunctionSpace>}
1459 jgs 108 """
1460 jgs 149 if self.reduceEquationOrder():
1461     return escript.ReducedSolution(self.getDomain())
1462     else:
1463     return escript.Solution(self.getDomain())
1464 jgs 108
1465     def getFunctionSpaceForSolution(self):
1466     """
1467 caltinay 2169 Returns the L{FunctionSpace<escript.FunctionSpace>} used to represent
1468     the solution.
1469 jgs 148
1470 jgs 149 @return: representation space of solution
1471     @rtype: L{FunctionSpace<escript.FunctionSpace>}
1472 jgs 108 """
1473 jgs 149 if self.reduceSolutionOrder():
1474     return escript.ReducedSolution(self.getDomain())
1475     else:
1476     return escript.Solution(self.getDomain())
1477 jgs 108
1478 caltinay 2169 # ==========================================================================
1479 jgs 148 # solver settings:
1480 caltinay 2169 # ==========================================================================
1481 gross 2474 def setSolverOptions(self,options=None):
1482 jgs 102 """
1483 gross 2474 Sets the solver options.
1484 jgs 148
1485 gross 2474 @param options: the new solver options. If equal C{None}, the solver options are set to the default.
1486     @type options: L{SolverOptions} or C{None}
1487     @note: The symmetry flag of options is overwritten by the symmetry flag of the L{LinearProblem}.
1488 jgs 102 """
1489 gross 2474 if options==None:
1490     self.__solver_options=SolverOptions()
1491     elif isinstance(options, SolverOptions):
1492     self.__solver_options=options
1493     else:
1494     raise ValueError,"options must be a SolverOptions object."
1495     self.__solver_options.setSymmetry(self.__sym)
1496    
1497     def getSolverOptions(self):
1498 jgs 148 """
1499 gross 2474 Returns the solver options
1500    
1501     @rtype: L{SolverOptions}
1502 jgs 148 """
1503 gross 2474 self.__solver_options.setSymmetry(self.__sym)
1504     return self.__solver_options
1505    
1506 jgs 148 def isUsingLumping(self):
1507     """
1508 caltinay 2169 Checks if matrix lumping is the current solver method.
1509 jgs 148
1510 caltinay 2169 @return: True if the current solver method is lumping
1511 jgs 148 @rtype: C{bool}
1512     """
1513 gross 2474 return self.getSolverOptions().getSolverMethod()==self.getSolverOptions().LUMPING
1514 caltinay 2169 # ==========================================================================
1515 jgs 148 # symmetry flag:
1516 caltinay 2169 # ==========================================================================
1517     def isSymmetric(self):
1518 jgs 102 """
1519 caltinay 2169 Checks if symmetry is indicated.
1520 jgs 149
1521 caltinay 2169 @return: True if a symmetric PDE is indicated, False otherwise
1522 jgs 149 @rtype: C{bool}
1523 gross 2474 @note: the method is equivalent to use getSolverOptions().isSymmetric()
1524 jgs 102 """
1525 gross 2474 self.getSolverOptions().isSymmetric()
1526 jgs 102
1527 caltinay 2169 def setSymmetryOn(self):
1528 jgs 102 """
1529 gross 2474 Sets the symmetry flag.
1530     @note: The method overwrites the symmetry flag set by the solver options
1531 jgs 102 """
1532 gross 2474 self.__sym=True
1533     self.getSolverOptions().setSymmetryOn()
1534 jgs 102
1535 caltinay 2169 def setSymmetryOff(self):
1536 jgs 102 """
1537 caltinay 2169 Clears the symmetry flag.
1538 gross 2474 @note: The method overwrites the symmetry flag set by the solver options
1539 jgs 102 """
1540 gross 2474 self.__sym=False
1541     self.getSolverOptions().setSymmetryOff()
1542 jgs 102
1543 gross 2474 def setSymmetry(self,flag=False):
1544 jgs 148 """
1545 caltinay 2169 Sets the symmetry flag to C{flag}.
1546 jgs 149
1547 caltinay 2169 @param flag: If True, the symmetry flag is set otherwise reset.
1548 jgs 148 @type flag: C{bool}
1549 gross 2474 @note: The method overwrites the symmetry flag set by the solver options
1550 jgs 148 """
1551 gross 2474 self.getSolverOptions().setSymmetry(flag)
1552 caltinay 2169 # ==========================================================================
1553 jgs 148 # function space handling for the equation as well as the solution
1554 caltinay 2169 # ==========================================================================
1555     def setReducedOrderOn(self):
1556 jgs 102 """
1557 caltinay 2169 Switches reduced order on for solution and equation representation.
1558 jgs 149
1559 caltinay 2169 @raise RuntimeError: if order reduction is altered after a coefficient has
1560     been set
1561 jgs 102 """
1562     self.setReducedOrderForSolutionOn()
1563     self.setReducedOrderForEquationOn()
1564    
1565 caltinay 2169 def setReducedOrderOff(self):
1566 jgs 102 """
1567 caltinay 2169 Switches reduced order off for solution and equation representation
1568 jgs 149
1569 caltinay 2169 @raise RuntimeError: if order reduction is altered after a coefficient has
1570     been set
1571 jgs 102 """
1572     self.setReducedOrderForSolutionOff()
1573     self.setReducedOrderForEquationOff()
1574    
1575 caltinay 2169 def setReducedOrderTo(self,flag=False):
1576 jgs 102 """
1577 caltinay 2169 Sets order reduction state for both solution and equation representation
1578     according to flag.
1579    
1580     @param flag: if True, the order reduction is switched on for both solution
1581     and equation representation, otherwise or if flag is not
1582     present order reduction is switched off
1583 jgs 148 @type flag: C{bool}
1584 caltinay 2169 @raise RuntimeError: if order reduction is altered after a coefficient has
1585     been set
1586 jgs 102 """
1587     self.setReducedOrderForSolutionTo(flag)
1588     self.setReducedOrderForEquationTo(flag)
1589    
1590 jgs 148
1591 caltinay 2169 def setReducedOrderForSolutionOn(self):
1592 jgs 102 """
1593 caltinay 2169 Switches reduced order on for solution representation.
1594 jgs 149
1595 caltinay 2169 @raise RuntimeError: if order reduction is altered after a coefficient has
1596     been set
1597 jgs 102 """
1598 jgs 149 if not self.__reduce_solution_order:
1599     if self.__altered_coefficients:
1600     raise RuntimeError,"order cannot be altered after coefficients have been defined."
1601 caltinay 2169 self.trace("Reduced order is used for solution representation.")
1602 jgs 149 self.__reduce_solution_order=True
1603 gross 1841 self.initializeSystem()
1604 jgs 102
1605 caltinay 2169 def setReducedOrderForSolutionOff(self):
1606 jgs 102 """
1607 caltinay 2169 Switches reduced order off for solution representation
1608 jgs 149
1609 caltinay 2169 @raise RuntimeError: if order reduction is altered after a coefficient has
1610     been set.
1611 jgs 102 """
1612 jgs 149 if self.__reduce_solution_order:
1613     if self.__altered_coefficients:
1614     raise RuntimeError,"order cannot be altered after coefficients have been defined."
1615 jgs 148 self.trace("Full order is used to interpolate solution.")
1616 jgs 149 self.__reduce_solution_order=False
1617 gross 1841 self.initializeSystem()
1618 jgs 102
1619 caltinay 2169 def setReducedOrderForSolutionTo(self,flag=False):
1620 jgs 102 """
1621 caltinay 2169 Sets order reduction state for solution representation according to flag.
1622 jgs 102
1623 caltinay 2169 @param flag: if flag is True, the order reduction is switched on for
1624     solution representation, otherwise or if flag is not present
1625     order reduction is switched off
1626 jgs 148 @type flag: C{bool}
1627 caltinay 2169 @raise RuntimeError: if order reduction is altered after a coefficient has
1628     been set
1629 jgs 102 """
1630     if flag:
1631     self.setReducedOrderForSolutionOn()
1632     else:
1633     self.setReducedOrderForSolutionOff()
1634 jgs 148
1635 caltinay 2169 def setReducedOrderForEquationOn(self):
1636 jgs 102 """
1637 caltinay 2169 Switches reduced order on for equation representation.
1638 jgs 149
1639 caltinay 2169 @raise RuntimeError: if order reduction is altered after a coefficient has
1640     been set
1641 jgs 102 """
1642 jgs 149 if not self.__reduce_equation_order:
1643     if self.__altered_coefficients:
1644     raise RuntimeError,"order cannot be altered after coefficients have been defined."
1645 jgs 148 self.trace("Reduced order is used for test functions.")
1646 jgs 149 self.__reduce_equation_order=True
1647 gross 1841 self.initializeSystem()
1648 jgs 102
1649 caltinay 2169 def setReducedOrderForEquationOff(self):
1650 jgs 102 """
1651 caltinay 2169 Switches reduced order off for equation representation.
1652 jgs 149
1653 caltinay 2169 @raise RuntimeError: if order reduction is altered after a coefficient has
1654     been set
1655 jgs 102 """
1656 jgs 149 if self.__reduce_equation_order:
1657     if self.__altered_coefficients:
1658     raise RuntimeError,"order cannot be altered after coefficients have been defined."
1659 jgs 148 self.trace("Full order is used for test functions.")
1660 jgs 149 self.__reduce_equation_order=False
1661 gross 1841 self.initializeSystem()
1662 jgs 102
1663 caltinay 2169 def setReducedOrderForEquationTo(self,flag=False):
1664 jgs 102 """
1665 caltinay 2169 Sets order reduction state for equation representation according to flag.
1666 jgs 102
1667 caltinay 2169 @param flag: if flag is True, the order reduction is switched on for
1668     equation representation, otherwise or if flag is not present
1669     order reduction is switched off
1670 jgs 148 @type flag: C{bool}
1671 caltinay 2169 @raise RuntimeError: if order reduction is altered after a coefficient has
1672     been set
1673 jgs 102 """
1674     if flag:
1675     self.setReducedOrderForEquationOn()
1676     else:
1677     self.setReducedOrderForEquationOff()
1678 gross 2474 def getOperatorType(self):
1679 gross 1841 """
1680 caltinay 2169 Returns the current system type.
1681 gross 1841 """
1682 gross 2474 return self.__operator_type
1683 jgs 148
1684 gross 1841 def checkSymmetricTensor(self,name,verbose=True):
1685     """
1686 caltinay 2169 Tests a coefficient for symmetry.
1687 jgs 148
1688 gross 1841 @param name: name of the coefficient
1689     @type name: C{str}
1690 caltinay 2169 @param verbose: if set to True or not present a report on coefficients
1691     which break the symmetry is printed.
1692 gross 1841 @type verbose: C{bool}
1693 caltinay 2169 @return: True if coefficient C{name} is symmetric
1694 gross 1841 @rtype: C{bool}
1695     """
1696 gross 2474 SMALL_TOLERANCE=util.EPSILON*10.
1697 gross 1841 A=self.getCoefficient(name)
1698     verbose=verbose or self.__debug
1699     out=True
1700     if not A.isEmpty():
1701 gross 2474 tol=util.Lsup(A)*SMALL_TOLERANCE
1702 gross 1841 s=A.getShape()
1703     if A.getRank() == 4:
1704     if s[0]==s[2] and s[1] == s[3]:
1705     for i in range(s[0]):
1706     for j in range(s[1]):
1707     for k in range(s[2]):
1708     for l in range(s[3]):
1709     if util.Lsup(A[i,j,k,l]-A[k,l,i,j])>tol:
1710     if verbose: print "non-symmetric problem as %s[%d,%d,%d,%d]!=%s[%d,%d,%d,%d]"%(name,i,j,k,l,name,k,l,i,j)
1711     out=False
1712     else:
1713     if verbose: print "non-symmetric problem because of inappropriate shape %s of coefficient %s."%(s,name)
1714     out=False
1715     elif A.getRank() == 2:
1716     if s[0]==s[1]:
1717     for j in range(s[0]):
1718     for l in range(s[1]):
1719     if util.Lsup(A[j,l]-A[l,j])>tol:
1720     if verbose: print "non-symmetric problem because %s[%d,%d]!=%s[%d,%d]"%(name,j,l,name,l,j)
1721     out=False
1722     else:
1723     if verbose: print "non-symmetric problem because of inappropriate shape %s of coefficient %s."%(s,name)
1724     out=False
1725     elif A.getRank() == 0:
1726     pass
1727     else:
1728     raise ValueError,"Cannot check rank %s of %s."%(A.getRank(),name)
1729     return out
1730 jgs 148
1731 gross 1841 def checkReciprocalSymmetry(self,name0,name1,verbose=True):
1732     """
1733 caltinay 2169 Tests two coefficients for reciprocal symmetry.
1734 jgs 148
1735 gross 1841 @param name0: name of the first coefficient
1736 jfenwick 2061 @type name0: C{str}
1737 gross 1841 @param name1: name of the second coefficient
1738 jfenwick 2061 @type name1: C{str}
1739 caltinay 2169 @param verbose: if set to True or not present a report on coefficients
1740     which break the symmetry is printed
1741 gross 1841 @type verbose: C{bool}
1742 caltinay 2169 @return: True if coefficients C{name0} and C{name1} are reciprocally
1743     symmetric.
1744 gross 1841 @rtype: C{bool}
1745     """
1746 gross 2474 SMALL_TOLERANCE=util.EPSILON*10.
1747 gross 1841 B=self.getCoefficient(name0)
1748     C=self.getCoefficient(name1)
1749     verbose=verbose or self.__debug
1750     out=True
1751     if B.isEmpty() and not C.isEmpty():
1752     if verbose: print "non-symmetric problem because %s is not present but %s is"%(name0,name1)
1753     out=False
1754     elif not B.isEmpty() and C.isEmpty():
1755     if verbose: print "non-symmetric problem because %s is not present but %s is"%(name0,name1)
1756     out=False
1757     elif not B.isEmpty() and not C.isEmpty():
1758     sB=B.getShape()
1759     sC=C.getShape()
1760 gross 2474 tol=(util.Lsup(B)+util.Lsup(C))*SMALL_TOLERANCE/2.
1761 gross 1841 if len(sB) != len(sC):
1762     if verbose: print "non-symmetric problem because ranks of %s (=%s) and %s (=%s) are different."%(name0,len(sB),name1,len(sC))
1763     out=False
1764     else:
1765     if len(sB)==0:
1766     if util.Lsup(B-C)>tol:
1767     if verbose: print "non-symmetric problem because %s!=%s"%(name0,name1)
1768     out=False
1769     elif len(sB)==1:
1770     if sB[0]==sC[0]:
1771     for j in range(sB[0]):
1772     if util.Lsup(B[j]-C[j])>tol:
1773     if verbose: print "non-symmetric PDE because %s[%d]!=%s[%d]"%(name0,j,name1,j)
1774     out=False
1775     else:
1776     if verbose: print "non-symmetric problem because of inappropriate shapes %s and %s of coefficients %s and %s, respectively."%(sB,sC,name0,name1)
1777     elif len(sB)==3:
1778     if sB[0]==sC[1] and sB[1]==sC[2] and sB[2]==sC[0]:
1779     for i in range(sB[0]):
1780     for j in range(sB[1]):
1781     for k in range(sB[2]):
1782     if util.Lsup(B[i,j,k]-C[k,i,j])>tol:
1783 gross 2100 if verbose: print "non-symmetric problem because %s[%d,%d,%d]!=%s[%d,%d,%d]"%(name0,i,j,k,name1,k,i,j)
1784 gross 1841 out=False
1785     else:
1786     if verbose: print "non-symmetric problem because of inappropriate shapes %s and %s of coefficients %s and %s, respectively."%(sB,sC,name0,name1)
1787 jgs 148 else:
1788 gross 1841 raise ValueError,"Cannot check rank %s of %s and %s."%(len(sB),name0,name1)
1789     return out
1790 jgs 148
1791 caltinay 2169 def getCoefficient(self,name):
1792 jgs 108 """
1793 caltinay 2169 Returns the value of the coefficient C{name}.
1794 jgs 108
1795 caltinay 2169 @param name: name of the coefficient requested
1796 jgs 148 @type name: C{string}
1797 caltinay 2169 @return: the value of the coefficient
1798 jgs 149 @rtype: L{Data<escript.Data>}
1799 caltinay 2169 @raise IllegalCoefficient: if C{name} is not a coefficient of the PDE
1800 jgs 148 """
1801     if self.hasCoefficient(name):
1802 gross 1841 return self.__COEFFICIENTS[name].getValue()
1803 jgs 148 else:
1804     raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1805 jgs 108
1806 caltinay 2169 def hasCoefficient(self,name):
1807 jgs 148 """
1808 caltinay 2169 Returns True if C{name} is the name of a coefficient.
1809 jgs 108
1810 caltinay 2169 @param name: name of the coefficient enquired
1811 jgs 148 @type name: C{string}
1812 caltinay 2169 @return: True if C{name} is the name of a coefficient of the general PDE,
1813     False otherwise
1814 jgs 149 @rtype: C{bool}
1815 jgs 148 """
1816 gross 1841 return self.__COEFFICIENTS.has_key(name)
1817 jgs 108
1818 caltinay 2169 def createCoefficient(self, name):
1819 jgs 148 """
1820 caltinay 2169 Creates a L{Data<escript.Data>} object corresponding to coefficient
1821     C{name}.
1822 jgs 108
1823 caltinay 2169 @return: the coefficient C{name} initialized to 0
1824 jgs 149 @rtype: L{Data<escript.Data>}
1825 caltinay 2169 @raise IllegalCoefficient: if C{name} is not a coefficient of the PDE
1826 jgs 148 """
1827     if self.hasCoefficient(name):
1828     return escript.Data(0.,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name))
1829     else:
1830     raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1831    
1832 caltinay 2169 def getFunctionSpaceForCoefficient(self,name):
1833 jgs 148 """
1834 caltinay 2169 Returns the L{FunctionSpace<escript.FunctionSpace>} to be used for
1835     coefficient C{name}.
1836 jgs 148
1837 caltinay 2169 @param name: name of the coefficient enquired
1838 jgs 148 @type name: C{string}
1839 caltinay 2169 @return: the function space to be used for coefficient C{name}
1840 jgs 149 @rtype: L{FunctionSpace<escript.FunctionSpace>}
1841 caltinay 2169 @raise IllegalCoefficient: if C{name} is not a coefficient of the PDE
1842 jgs 148 """
1843     if self.hasCoefficient(name):
1844 gross 1841 return self.__COEFFICIENTS[name].getFunctionSpace(self.getDomain())
1845 jgs 148 else:
1846     raise ValueError,"unknown coefficient %s requested"%name
1847 gross 1841
1848 caltinay 2169 def getShapeOfCoefficient(self,name):
1849 jgs 148 """
1850 caltinay 2169 Returns the shape of the coefficient C{name}.
1851 jgs 148
1852 caltinay 2169 @param name: name of the coefficient enquired
1853 jgs 148 @type name: C{string}
1854 caltinay 2169 @return: the shape of the coefficient C{name}
1855 jgs 149 @rtype: C{tuple} of C{int}
1856 caltinay 2169 @raise IllegalCoefficient: if C{name} is not a coefficient of the PDE
1857 jgs 148 """
1858     if self.hasCoefficient(name):
1859 gross 1841 return self.__COEFFICIENTS[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())
1860 jgs 148 else:
1861     raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1862    
1863 caltinay 2169 def resetAllCoefficients(self):
1864 jgs 148 """
1865 caltinay 2169 Resets all coefficients to their default values.
1866 jgs 148 """
1867 gross 1841 for i in self.__COEFFICIENTS.iterkeys():
1868     self.__COEFFICIENTS[i].resetValue()
1869 jgs 148
1870 caltinay 2169 def alteredCoefficient(self,name):
1871 jgs 148 """
1872 caltinay 2169 Announces that coefficient C{name} has been changed.
1873 jgs 148
1874 caltinay 2169 @param name: name of the coefficient affected
1875 jgs 148 @type name: C{string}
1876 caltinay 2169 @raise IllegalCoefficient: if C{name} is not a coefficient of the PDE
1877     @note: if C{name} is q or r, the method will not trigger a rebuild of the
1878     system as constraints are applied to the solved system.
1879 jgs 148 """
1880     if self.hasCoefficient(name):
1881     self.trace("Coefficient %s has been altered."%name)
1882 jgs 149 if not ((name=="q" or name=="r") and self.isUsingLumping()):
1883 gross 1841 if self.__COEFFICIENTS[name].isAlteringOperator(): self.invalidateOperator()
1884     if self.__COEFFICIENTS[name].isAlteringRightHandSide(): self.invalidateRightHandSide()
1885 jgs 148 else:
1886     raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1887    
1888 caltinay 2169 def validSolution(self):
1889 gross 1841 """
1890 caltinay 2169 Marks the solution as valid.
1891 gross 1841 """
1892     self.__is_solution_valid=True
1893 jgs 149
1894 caltinay 2169 def invalidateSolution(self):
1895 gross 1841 """
1896 caltinay 2169 Indicates the PDE has to be resolved if the solution is requested.
1897 gross 1841 """
1898     self.trace("System will be resolved.")
1899     self.__is_solution_valid=False
1900 jgs 148
1901 gross 1841 def isSolutionValid(self):
1902     """
1903 caltinay 2169 Returns True if the solution is still valid.
1904 gross 1841 """
1905 gross 2535 if not self.getDomainStatus()==self.getSystemStatus(): self.invalidateSolution()
1906 gross 2474 if self.__solution_rtol>self.getSolverOptions().getTolerance() or \
1907     self.__solution_atol>self.getSolverOptions().getAbsoluteTolerance():
1908     self.invalidateSolution()
1909 gross 1841 return self.__is_solution_valid
1910    
1911     def validOperator(self):
1912     """
1913 caltinay 2169 Marks the operator as valid.
1914 gross 1841 """
1915     self.__is_operator_valid=True
1916    
1917 caltinay 2169 def invalidateOperator(self):
1918 gross 1841 """
1919 caltinay 2169 Indicates the operator has to be rebuilt next time it is used.
1920 gross 1841 """
1921     self.trace("Operator will be rebuilt.")
1922     self.invalidateSolution()
1923     self.__is_operator_valid=False
1924    
1925     def isOperatorValid(self):
1926     """
1927 caltinay 2169 Returns True if the operator is still valid.
1928 gross 1841 """
1929 gross 2535 if not self.getDomainStatus()==self.getSystemStatus(): self.invalidateOperator()
1930     if not self.getRequiredOperatorType()==self.getOperatorType(): self.invalidateOperator()
1931 gross 1841 return self.__is_operator_valid
1932    
1933     def validRightHandSide(self):
1934     """
1935 caltinay 2169 Marks the right hand side as valid.
1936 gross 1841 """
1937     self.__is_RHS_valid=True
1938    
1939 caltinay 2169 def invalidateRightHandSide(self):
1940 gross 1841 """
1941 caltinay 2169 Indicates the right hand side has to be rebuilt next time it is used.
1942 gross 1841 """
1943 gross 2535 self.trace("Right hand side has to be rebuilt.")
1944 gross 1841 self.invalidateSolution()
1945     self.__is_RHS_valid=False
1946    
1947     def isRightHandSideValid(self):
1948     """
1949 caltinay 2169 Returns True if the operator is still valid.
1950 gross 1841 """
1951 gross 2535 if not self.getDomainStatus()==self.getSystemStatus(): self.invalidateRightHandSide()
1952 gross 1841 return self.__is_RHS_valid
1953    
1954 caltinay 2169 def invalidateSystem(self):
1955 gross 1841 """
1956 caltinay 2169 Announces that everything has to be rebuilt.
1957 gross 1841 """
1958     self.invalidateSolution()
1959     self.invalidateOperator()
1960     self.invalidateRightHandSide()
1961    
1962     def isSystemValid(self):
1963     """
1964 caltinay 2169 Returns True if the system (including solution) is still vaild.
1965 gross 1841 """
1966     return self.isSolutionValid() and self.isOperatorValid() and self.isRightHandSideValid()
1967    
1968 caltinay 2169 def initializeSystem(self):
1969 gross 1841 """
1970 caltinay 2169 Resets the system clearing the operator, right hand side and solution.
1971 gross 1841 """
1972     self.trace("New System has been created.")
1973 gross 2474 self.__operator_type=None
1974 gross 2535 self.setSystemStatus()
1975 gross 1841 self.__operator=escript.Operator()
1976     self.__righthandside=escript.Data()
1977     self.__solution=escript.Data()
1978     self.invalidateSystem()
1979    
1980 caltinay 2169 def getOperator(self):
1981 gross 1841 """
1982 caltinay 2169 Returns the operator of the linear problem.
1983 gross 1841
1984     @return: the operator of the problem
1985     """
1986     return self.getSystem()[0]
1987    
1988 caltinay 2169 def getRightHandSide(self):
1989 gross 1841 """
1990 caltinay 2169 Returns the right hand side of the linear problem.
1991 gross 1841
1992     @return: the right hand side of the problem
1993     @rtype: L{Data<escript.Data>}
1994     """
1995     return self.getSystem()[1]
1996    
1997 caltinay 2169 def createRightHandSide(self):
1998 gross 1841 """
1999 caltinay 2169 Returns an instance of a new right hand side.
2000 gross 1841 """
2001     self.trace("New right hand side is allocated.")
2002     if self.getNumEquations()>1:
2003     return escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True)
2004     else:
2005     return escript.Data(0.,(),self.getFunctionSpaceForEquation(),True)
2006    
2007 caltinay 2169 def createSolution(self):
2008 gross 1841 """
2009 caltinay 2169 Returns an instance of a new solution.
2010 gross 1841 """
2011     self.trace("New solution is allocated.")
2012     if self.getNumSolutions()>1:
2013     return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)
2014     else:
2015     return escript.Data(0.,(),self.getFunctionSpaceForSolution(),True)
2016    
2017 caltinay 2169 def resetSolution(self):
2018 gross 1841 """
2019 caltinay 2169 Sets the solution to zero.
2020 gross 1841 """
2021     if self.__solution.isEmpty():
2022     self.__solution=self.createSolution()
2023     else:
2024     self.__solution.setToZero()
2025     self.trace("Solution is reset to zero.")
2026 caltinay 2169
2027 gross 1841 def setSolution(self,u):
2028     """
2029 gross 2474 Sets the solution assuming that makes the system valid with the tolrance
2030     defined by the solver options
2031 gross 1841 """
2032 gross 2474 self.__solution_rtol=self.getSolverOptions().getTolerance()
2033     self.__solution_atol=self.getSolverOptions().getAbsoluteTolerance()
2034 gross 1841 self.__solution=u
2035     self.validSolution()
2036    
2037     def getCurrentSolution(self):
2038     """
2039 caltinay 2169 Returns the solution in its current state.
2040 gross 1841 """
2041     if self.__solution.isEmpty(): self.__solution=self.createSolution()
2042     return self.__solution
2043    
2044 caltinay 2169 def resetRightHandSide(self):
2045 gross 1841 """
2046 caltinay 2169 Sets the right hand side to zero.
2047 gross 1841 """
2048     if self.__righthandside.isEmpty():
2049     self.__righthandside=self.createRightHandSide()
2050     else:
2051     self.__righthandside.setToZero()
2052     self.trace("Right hand side is reset to zero.")
2053    
2054     def getCurrentRightHandSide(self):
2055     """
2056 caltinay 2169 Returns the right hand side in its current state.
2057 gross 1841 """
2058     if self.__righthandside.isEmpty(): self.__righthandside=self.createRightHandSide()
2059     return self.__righthandside
2060    
2061 caltinay 2169 def resetOperator(self):
2062 gross 1841 """
2063 caltinay 2169 Makes sure that the operator is instantiated and returns it initialized
2064     with zeros.
2065 gross 1841 """
2066 gross 2474 if self.getOperatorType() == None:
2067 gross 1841 if self.isUsingLumping():
2068     self.__operator=self.createSolution()
2069     else:
2070     self.__operator=self.createOperator()
2071 gross 2474 self.__operator_type=self.getRequiredOperatorType()
2072 gross 1841 else:
2073     if self.isUsingLumping():
2074     self.__operator.setToZero()
2075     else:
2076 gross 2551 if self.getOperatorType() == self.getRequiredOperatorType():
2077     self.__operator.resetValues()
2078     else:
2079     self.__operator=self.createOperator()
2080     self.__operator_type=self.getRequiredOperatorType()
2081 gross 1841 self.trace("Operator reset to zero")
2082    
2083     def getCurrentOperator(self):
2084     """
2085 caltinay 2169 Returns the operator in its current state.
2086 gross 1841 """
2087     return self.__operator
2088    
2089 caltinay 2169 def setValue(self,**coefficients):
2090 jgs 148 """
2091 caltinay 2169 Sets new values to coefficients.
2092 jgs 148
2093 caltinay 2169 @raise IllegalCoefficient: if an unknown coefficient keyword is used
2094 jgs 148 """
2095 jgs 108 # check if the coefficients are legal:
2096     for i in coefficients.iterkeys():
2097     if not self.hasCoefficient(i):
2098 jgs 148 raise IllegalCoefficient,"Attempt to set unknown coefficient %s"%i
2099 jgs 108 # if the number of unknowns or equations is still unknown we try to estimate them:
2100 jgs 148 if self.__numEquations==None or self.__numSolutions==None:
2101 jgs 108 for i,d in coefficients.iteritems():
2102     if hasattr(d,"shape"):
2103     s=d.shape
2104     elif hasattr(d,"getShape"):
2105     s=d.getShape()
2106     else:
2107 jfenwick 2455 s=numpy.array(d).shape
2108 jgs 108 if s!=None:
2109     # get number of equations and number of unknowns:
2110 gross 1841 res=self.__COEFFICIENTS[i].estimateNumEquationsAndNumSolutions(self.getDomain(),s)
2111 jgs 108 if res==None:
2112 jgs 148 raise IllegalCoefficientValue,"Illegal shape %s of coefficient %s"%(s,i)
2113 jgs 108 else:
2114 jgs 148 if self.__numEquations==None: self.__numEquations=res[0]
2115     if self.__numSolutions==None: self.__numSolutions=res[1]
2116 caltinay 2169 if self.__numEquations==None: raise UndefinedPDEError,"unidentified number of equations"
2117     if self.__numSolutions==None: raise UndefinedPDEError,"unidentified number of solutions"
2118 jgs 108 # now we check the shape of the coefficient if numEquations and numSolutions are set:
2119     for i,d in coefficients.iteritems():
2120 jgs 148 try:
2121 gross 1841 self.__COEFFICIENTS[i].setValue(self.getDomain(),
2122 caltinay 2169 self.getNumEquations(),self.getNumSolutions(),
2123     self.reduceEquationOrder(),self.reduceSolutionOrder(),d)
2124 gross 1072 self.alteredCoefficient(i)
2125     except IllegalCoefficientFunctionSpace,m:
2126     # if the function space is wrong then we try the reduced version:
2127     i_red=i+"_reduced"
2128 gross 1841 if (not i_red in coefficients.keys()) and i_red in self.__COEFFICIENTS.keys():
2129 gross 1072 try:
2130 gross 1841 self.__COEFFICIENTS[i_red].setValue(self.getDomain(),
2131 gross 1072 self.getNumEquations(),self.getNumSolutions(),
2132     self.reduceEquationOrder(),self.reduceSolutionOrder(),d)
2133     self.alteredCoefficient(i_red)
2134     except IllegalCoefficientValue,m:
2135     raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
2136     except IllegalCoefficientFunctionSpace,m:
2137     raise IllegalCoefficientFunctionSpace("Coefficient %s:%s"%(i,m))
2138     else:
2139     raise IllegalCoefficientFunctionSpace("Coefficient %s:%s"%(i,m))
2140 jgs 148 except IllegalCoefficientValue,m:
2141     raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
2142 jgs 149 self.__altered_coefficients=True
2143 jgs 108
2144 caltinay 2169 # ==========================================================================
2145     # methods that are typically overwritten when implementing a particular
2146     # linear problem
2147     # ==========================================================================
2148 gross 2474 def getRequiredOperatorType(self):
2149 gross 1841 """
2150 caltinay 2169 Returns the system type which needs to be used by the current set up.
2151 gross 1841
2152 caltinay 2169 @note: Typically this method is overwritten when implementing a
2153     particular linear problem.
2154 gross 1841 """
2155 gross 2474 return None
2156 gross 1841
2157 caltinay 2169 def createOperator(self):
2158 jgs 108 """
2159 caltinay 2169 Returns an instance of a new operator.
2160 gross 1841
2161 caltinay 2169 @note: This method is overwritten when implementing a particular
2162     linear problem.
2163 gross 1841 """
2164     return escript.Operator()
2165    
2166     def checkSymmetry(self,verbose=True):
2167     """
2168 caltinay 2169 Tests the PDE for symmetry.
2169 gross 1841
2170 caltinay 2169 @param verbose: if set to True or not present a report on coefficients
2171     which break the symmetry is printed
2172 gross 1841 @type verbose: C{bool}
2173 caltinay 2169 @return: True if the problem is symmetric
2174 gross 1841 @rtype: C{bool}
2175 caltinay 2169 @note: Typically this method is overwritten when implementing a
2176     particular linear problem.
2177 gross 1841 """
2178     out=True
2179     return out
2180    
2181 caltinay 2169 def getSolution(self,**options):
2182 gross 1841 """
2183 caltinay 2169 Returns the solution of the problem.
2184    
2185 gross 1841 @return: the solution
2186     @rtype: L{Data<escript.Data>}
2187    
2188 caltinay 2169 @note: This method is overwritten when implementing a particular
2189     linear problem.
2190 gross 1841 """
2191     return self.getCurrentSolution()
2192    
2193 caltinay 2169 def getSystem(self):
2194 gross 1841 """
2195 caltinay 2169 Returns the operator and right hand side of the PDE.
2196 jgs 149
2197     @return: the discrete version of the PDE
2198     @rtype: C{tuple} of L{Operator,<escript.Operator>} and L{Data<escript.Data>}.
2199 gross 1841
2200 caltinay 2169 @note: This method is overwritten when implementing a particular
2201     linear problem.
2202 jgs 108 """
2203 gross 1841 return (self.getCurrentOperator(), self.getCurrentRightHandSide())
2204    
2205     class LinearPDE(LinearProblem):
2206     """
2207     This class is used to define a general linear, steady, second order PDE
2208 caltinay 2169 for an unknown function M{u} on a given domain defined through a
2209     L{Domain<escript.Domain>} object.
2210 gross 1841
2211 caltinay 2169 For a single PDE having a solution with a single component the linear PDE
2212     is defined in the following form:
2213 gross 1841
2214     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)}
2215    
2216 caltinay 2169 where M{grad(F)} denotes the spatial derivative of M{F}. Einstein's
2217     summation convention, ie. summation over indexes appearing twice in a term
2218     of a sum performed, is used.
2219     The coefficients M{A}, M{B}, M{C}, M{D}, M{X} and M{Y} have to be specified
2220     through L{Data<escript.Data>} objects in L{Function<escript.Function>} and
2221     the coefficients M{A_reduced}, M{B_reduced}, M{C_reduced}, M{D_reduced},
2222     M{X_reduced} and M{Y_reduced} have to be specified through
2223     L{Data<escript.Data>} objects in L{ReducedFunction<escript.ReducedFunction>}.
2224     It is also allowed to use objects that can be converted into such
2225     L{Data<escript.Data>} objects. M{A} and M{A_reduced} are rank two, M{B},
2226     M{C}, M{X}, M{B_reduced}, M{C_reduced} and M{X_reduced} are rank one and
2227     M{D}, M{D_reduced}, M{Y} and M{Y_reduced} are scalar.
2228 gross 1841
2229     The following natural boundary conditions are considered:
2230    
2231     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}
2232    
2233 caltinay 2169 where M{n} is the outer normal field. Notice that the coefficients M{A},
2234     M{A_reduced}, M{B}, M{B_reduced}, M{X} and M{X_reduced} are defined in the
2235     PDE. The coefficients M{d} and M{y} are each a scalar in
2236     L{FunctionOnBoundary<escript.FunctionOnBoundary>} and the coefficients
2237     M{d_reduced} and M{y_reduced} are each a scalar in
2238     L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
2239 gross 1841
2240 caltinay 2169 Constraints for the solution prescribe the value of the solution at certain
2241     locations in the domain. They have the form
2242 gross 1841
2243 caltinay 2169 M{u=r} where M{q>0}
2244 gross 1841
2245 caltinay 2169 M{r} and M{q} are each scalar where M{q} is the characteristic function
2246     defining where the constraint is applied. The constraints override any
2247     other condition set by the PDE or the boundary condition.
2248 gross 1841
2249     The PDE is symmetrical if
2250    
2251 caltinay 2169 M{A[i,j]=A[j,i]} and M{B[j]=C[j]} and M{A_reduced[i,j]=A_reduced[j,i]}
2252     and M{B_reduced[j]=C_reduced[j]}
2253 gross 1841
2254 caltinay 2169 For a system of PDEs and a solution with several components the PDE has the
2255     form
2256 gross 1841
2257     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] }
2258    
2259 caltinay 2169 M{A} and M{A_reduced} are of rank four, M{B}, M{B_reduced}, M{C} and
2260     M{C_reduced} are each of rank three, M{D}, M{D_reduced}, M{X_reduced} and
2261     M{X} are each of rank two and M{Y} and M{Y_reduced} are of rank one.
2262 gross 1841 The natural boundary conditions take the form:
2263    
2264     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]}
2265    
2266 caltinay 2169 The coefficient M{d} is of rank two and M{y} is of rank one both in
2267     L{FunctionOnBoundary<escript.FunctionOnBoundary>}. The coefficients
2268     M{d_reduced} is of rank two and M{y_reduced} is of rank one both in
2269     L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
2270 gross 1841
2271 caltinay 2169 Constraints take the form
2272 gross 1841
2273     M{u[i]=r[i]} where M{q[i]>0}
2274    
2275 caltinay 2169 M{r} and M{q} are each rank one. Notice that at some locations not
2276     necessarily all components must have a constraint.
2277 gross 1841
2278     The system of PDEs is symmetrical if
2279    
2280 caltinay 2169 - M{A[i,j,k,l]=A[k,l,i,j]}
2281     - M{A_reduced[i,j,k,l]=A_reduced[k,l,i,j]}
2282     - M{B[i,j,k]=C[k,i,j]}
2283     - M{B_reduced[i,j,k]=C_reduced[k,i,j]}
2284     - M{D[i,k]=D[i,k]}
2285     - M{D_reduced[i,k]=D_reduced[i,k]}
2286     - M{d[i,k]=d[k,i]}
2287     - M{d_reduced[i,k]=d_reduced[k,i]}
2288 gross 1841
2289 caltinay 2169 L{LinearPDE} also supports solution discontinuities over a contact region
2290     in the domain. To specify the conditions across the discontinuity we are
2291     using the generalised flux M{J} which, in the case of a system of PDEs
2292     and several components of the solution, is defined as
2293 gross 1841
2294     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]}
2295    
2296 caltinay 2169 For the case of single solution component and single PDE M{J} is defined as
2297 gross 1841
2298     M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[j]+(B[i]+B_reduced[i])*u-X[i]-X_reduced[i]}
2299    
2300 caltinay 2169 In the context of discontinuities M{n} denotes the normal on the
2301     discontinuity pointing from side 0 towards side 1 calculated from
2302     L{getNormal<escript.FunctionSpace.getNormal>} of L{FunctionOnContactZero<escript.FunctionOnContactZero>}.
2303     For a system of PDEs the contact condition takes the form
2304 gross 1841
2305     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]}
2306    
2307 caltinay 2169 where M{J0} and M{J1} are the fluxes on side 0 and side 1 of the
2308     discontinuity, respectively. M{jump(u)}, which is the difference of the
2309     solution at side 1 and at side 0, denotes the jump of M{u} across
2310     discontinuity along the normal calculated by L{jump<util.jump>}.
2311     The coefficient M{d_contact} is of rank two and M{y_contact} is of rank one
2312     both in L{FunctionOnContactZero<escript.FunctionOnContactZero>} or
2313     L{FunctionOnContactOne<escript.FunctionOnContactOne>}.
2314     The coefficient M{d_contact_reduced} is of rank two and M{y_contact_reduced}
2315     is of rank one both in L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>}
2316     or L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}.
2317     In case of a single PDE and a single component solution the contact
2318     condition takes the form
2319 gross 1841
2320     M{n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)}
2321    
2322 caltinay 2169 In this case the coefficient M{d_contact} and M{y_contact} are each scalar
2323     both in L{FunctionOnContactZero<escript.FunctionOnContactZero>} or
2324     L{FunctionOnContactOne<escript.FunctionOnContactOne>} and the coefficient
2325     M{d_contact_reduced} and M{y_contact_reduced} are each scalar both in
2326     L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or
2327     L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}.
2328 gross 1841
2329 caltinay 2169 Typical usage::
2330 gross 1841
2331 caltinay 2169 p = LinearPDE(dom)
2332     p.setValue(A=kronecker(dom), D=1, Y=0.5)
2333     u = p.getSolution()
2334 gross 1841
2335     """
2336    
2337     def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):
2338     """
2339 caltinay 2169 Initializes a new linear PDE.
2340 gross 1841
2341     @param domain: domain of the PDE
2342     @type domain: L{Domain<escript.Domain>}
2343 caltinay 2169 @param numEquations: number of equations. If C{None} the number of
2344     equations is extracted from the PDE coefficients.
2345     @param numSolutions: number of solution components. If C{None} the number
2346     of solution components is extracted from the PDE
2347     coefficients.
2348     @param debug: if True debug information is printed
2349 gross 1841
2350     """
2351     super(LinearPDE, self).__init__(domain,numEquations,numSolutions,debug)
2352     #
2353     # the coefficients of the PDE:
2354     #
2355     self.introduceCoefficients(
2356     A=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2357     B=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2358     C=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2359     D=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2360     X=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
2361     Y=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2362     d=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2363     y=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2364     d_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2365     y_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2366     A_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2367     B_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2368     C_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2369     D_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2370     X_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
2371     Y_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2372     d_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2373     y_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2374     d_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2375     y_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2376     r=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.RIGHTHANDSIDE),
2377     q=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.BOTH) )
2378    
2379     def __str__(self):
2380     """
2381 caltinay 2169 Returns the string representation of the PDE.
2382 gross 1841
2383     @return: a simple representation of the PDE
2384     @rtype: C{str}
2385     """
2386     return "<LinearPDE %d>"%id(self)
2387    
2388 gross 2474 def getRequiredOperatorType(self):
2389 gross 1841 """
2390 caltinay 2169 Returns the system type which needs to be used by the current set up.
2391 gross 1841 """
2392 gross 2474 solver_options=self.getSolverOptions()
2393     return self.getDomain().getSystemMatrixTypeId(solver_options.getSolverMethod(), solver_options.getPreconditioner(),solver_options.getPackage(), solver_options.isSymmetric())
2394 gross 1841
2395     def checkSymmetry(self,verbose=True):
2396     """
2397 caltinay 2169 Tests the PDE for symmetry.
2398 gross 1841
2399 caltinay 2169 @param verbose: if set to True or not present a report on coefficients
2400     which break the symmetry is printed.
2401 gross 1841 @type verbose: C{bool}
2402 caltinay 2169 @return: True if the PDE is symmetric
2403 gross 1841 @rtype: L{bool}
2404 caltinay 2169 @note: This is a very expensive operation. It should be used for
2405     degugging only! The symmetry flag is not altered.
2406 gross 1841 """
2407     out=True
2408     out=out and self.checkSymmetricTensor("A", verbose)
2409     out=out and self.checkSymmetricTensor("A_reduced", verbose)
2410     out=out and self.checkReciprocalSymmetry("B","C", verbose)
2411     out=out and self.checkReciprocalSymmetry("B_reduced","C_reduced", verbose)
2412     out=out and self.checkSymmetricTensor("D", verbose)
2413     out=out and self.checkSymmetricTensor("D_reduced", verbose)
2414     out=out and self.checkSymmetricTensor("d", verbose)
2415     out=out and self.checkSymmetricTensor("d_reduced", verbose)
2416     out=out and self.checkSymmetricTensor("d_contact", verbose)
2417     out=out and self.checkSymmetricTensor("d_contact_reduced", verbose)
2418     return out
2419    
2420 caltinay 2169 def createOperator(self):
2421 gross 1841 """
2422 caltinay 2169 Returns an instance of a new operator.
2423 gross 1841 """
2424 gross 2474 optype=self.getRequiredOperatorType()
2425     self.trace("New operator of type %s is allocated."%optype)
2426 gross 1841 return self.getDomain().newOperator( \
2427     self.getNumEquations(), \
2428     self.getFunctionSpaceForEquation(), \
2429     self.getNumSolutions(), \
2430     self.getFunctionSpaceForSolution(), \
2431 gross 2474 optype)
2432 gross 1841
2433 gross 2474 def getSolution(self):
2434 gross 1841 """
2435 caltinay 2169 Returns the solution of the PDE.
2436 gross 1841
2437     @return: the solution
2438     @rtype: L{Data<escript.Data>}
2439     """
2440 gross 2474 option_class=self.getSolverOptions()
2441 gross 1841 if not self.isSolutionValid():
2442     mat,f=self.getSystem()
2443 jgs 108 if self.isUsingLumping():
2444 gross 1841 self.setSolution(f*1/mat)
2445     else:
2446     self.trace("PDE is resolved.")
2447 gross 2474 self.trace("solver options: %s"%str(option_class))
2448     self.setSolution(mat.solve(f,option_class))
2449 gross 1841 return self.getCurrentSolution()
2450    
2451 caltinay 2169 def getSystem(self):
2452 gross 1841 """
2453 caltinay 2169 Returns the operator and right hand side of the PDE.
2454 gross 1841
2455     @return: the discrete version of the PDE
2456 caltinay 2169 @rtype: C{tuple} of L{Operator,<escript.Operator>} and
2457     L{Data<escript.Data>}
2458 gross 1841 """
2459     if not self.isOperatorValid() or not self.isRightHandSideValid():
2460     if self.isUsingLumping():
2461     if not self.isOperatorValid():
2462 caltinay 2169 if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():
2463 gross 531 raise TypeError,"Lumped matrix requires same order for equations and unknowns"
2464 caltinay 2169 if not self.getCoefficient("A").isEmpty():
2465 gross 531 raise ValueError,"coefficient A in lumped matrix may not be present."
2466 gross 1841 if not self.getCoefficient("B").isEmpty():
2467 gross 1072 raise ValueError,"coefficient B in lumped matrix may not be present."