/[escript]/branches/symbolic_from_3470/escript/py_src/linearPDEs.py
ViewVC logotype

Annotation of /branches/symbolic_from_3470/escript/py_src/linearPDEs.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3402 - (hide annotations)
Tue Dec 7 07:36:12 2010 UTC (8 years, 11 months ago) by gross
Original Path: trunk/escript/py_src/linearPDEs.py
File MIME type: text/x-python
File size: 161021 byte(s)
AMG support now block size > 3 (requires clapack or MKL).

additional diagnostics for AMG 


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