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