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