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