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