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