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