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