/[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 3434 - (hide annotations)
Fri Jan 7 03:35:58 2011 UTC (8 years, 10 months ago) by gross
Original Path: trunk/escript/py_src/linearPDEs.py
File MIME type: text/x-python
File size: 160482 byte(s)
initial value changed
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 gross 3433 def setMinCoarseMatrixSize(self,size=None):
415 artak 2524 """
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 3434 def setNumSweeps(self,sweeps=2):
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 gross 3434 def setNumPreSweeps(self,sweeps=2):
711 gross 2470 """
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 gross 3431 def setNumPostSweeps(self,sweeps=1):
729 gross 2470 """
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 sparsity<0:
1033     raise ValueError,"sparsity must be non-negative."
1034     if sparsity>1:
1035     raise ValueError,"sparsity must be less than 1."
1036     self.__min_sparsity=sparsity
1037 gross 3094
1038 gross 3283 def getMinCoarseMatrixSparsity(self):
1039     """
1040     Returns the minimum sparsity on the coarsest level. Typically
1041     a direct solver is used when the sparsity becomes bigger than
1042     the set limit.
1043    
1044     :return: minimual sparsity
1045     :rtype: ``float``
1046     """
1047     return self.__min_sparsity
1048    
1049     def setNumRefinements(self,refinements=2):
1050     """
1051     Sets the number of refinement steps to refine the solution when a direct solver is applied.
1052    
1053     :param refinements: number of refinements
1054     :type refinements: non-negative ``int``
1055     """
1056     refinements=int(refinements)
1057     if refinements<0:
1058     raise ValueError,"number of refinements must be non-negative."
1059     self.__refinements=refinements
1060    
1061     def getNumRefinements(self):
1062     """
1063     Returns the number of refinement steps to refine the solution when a direct solver is applied.
1064    
1065     :rtype: non-negative ``int``
1066     """
1067     return self.__refinements
1068    
1069     def setNumCoarseMatrixRefinements(self,refinements=2):
1070     """
1071     Sets the number of refinement steps to refine the solution on the coarset level when a direct solver is applied.
1072    
1073     :param refinements: number of refinements
1074     :type refinements: non-negative ``int``
1075     """
1076     refinements=int(refinements)
1077     if refinements<0:
1078     raise ValueError,"number of refinements must be non-negative."
1079     self.__coarse_refinements=refinements
1080    
1081     def getNumCoarseMatrixRefinements(self):
1082     """
1083     Returns the number of resfinement steps to refine the solution on the coarset level when a direct solver is applied.
1084    
1085     :rtype: non-negative ``int``
1086     """
1087     return self.__coarse_refinements
1088 gross 3352
1089     def usePanel(self):
1090     """
1091     Returns ``True`` if a panel is used to serach for unknown in the AMG coarsening, The panel approach is normally faster
1092     but can lead to larger coarse level systems.
1093    
1094     :return: ``True`` if a panel is used to find unknowns in AMG coarsening
1095     :rtype: ``bool``
1096     """
1097     return self.__use_panel
1098    
1099     def setUsePanelOn(self):
1100     """
1101     Sets the flag to use a panel to find unknowns in AMG coarsening
1102     """
1103     self.__use_panel=True
1104    
1105     def setUsePanelOff(self):
1106     """
1107     Sets the flag to use a panel to find unknowns in AMG coarsening to off
1108     """
1109     self.__use_panel=False
1110    
1111     def setUsePanel(self,use=True):
1112     """
1113     Sets the flag to use a panel to find unknowns in AMG coarsening
1114    
1115 jfenwick 3360 :param use: If ``True``,a panel is used to find unknowns in AMG coarsening
1116 gross 3352 :type use: ``bool``
1117     """
1118     if use:
1119     self.setUsePanelOn()
1120     else:
1121     self.setUsePanelOff()
1122    
1123     def useDirectInterpolation(self):
1124     """
1125     Returns ``True`` if direct interpolation is used in AMG.
1126    
1127     :return: ``True`` if direct interpolation is used in AMG.
1128     :rtype: ``bool``
1129     """
1130     return self.__use_direct_interpolation
1131    
1132     def setUseDirectInterpolationOn(self):
1133     """
1134     Sets the flag to use direct interpolation in AMG
1135     """
1136     self.__use_direct_interpolation=True
1137    
1138     def setUseDirectInterpolationOff(self):
1139     """
1140     Sets the flag to use direct interpolation in AMG
1141     """
1142     self.__use_direct_interpolation=False
1143    
1144     def setUseDirectInterpolation(self,use=False):
1145     """
1146     Sets the flag to use direct interpolation in AMG
1147    
1148     :param use: If ``True``, direct interpolation in AMG
1149     :type use: ``bool``
1150     """
1151     if use:
1152     self.setUseDirectInterpolationOn()
1153     else:
1154     self.setUseDirectInterpolationOff()
1155    
1156    
1157 jgs 148 class IllegalCoefficient(ValueError):
1158 jgs 102 """
1159 caltinay 2169 Exception that is raised if an illegal coefficient of the general or
1160     particular PDE is requested.
1161 jgs 148 """
1162 gross 1072 pass
1163 jgs 102
1164 jgs 148 class IllegalCoefficientValue(ValueError):
1165 jgs 102 """
1166 caltinay 2169 Exception that is raised if an incorrect value for a coefficient is used.
1167 jgs 148 """
1168 gross 1072 pass
1169 jgs 122
1170 gross 1072 class IllegalCoefficientFunctionSpace(ValueError):
1171     """
1172 caltinay 2169 Exception that is raised if an incorrect function space for a coefficient
1173     is used.
1174 gross 1072 """
1175    
1176 jgs 148 class UndefinedPDEError(ValueError):
1177     """
1178 caltinay 2169 Exception that is raised if a PDE is not fully defined yet.
1179 jgs 148 """
1180 gross 1072 pass
1181 jgs 102
1182 gross 1841 class PDECoef(object):
1183 jgs 102 """
1184 caltinay 2169 A class for describing a PDE coefficient.
1185 jgs 149
1186 jfenwick 2625 :cvar INTERIOR: indicator that coefficient is defined on the interior of
1187 caltinay 2169 the PDE domain
1188 jfenwick 2625 :cvar BOUNDARY: indicator that coefficient is defined on the boundary of
1189 caltinay 2169 the PDE domain
1190 jfenwick 2625 :cvar CONTACT: indicator that coefficient is defined on the contact region
1191 caltinay 2169 within the PDE domain
1192 jfenwick 2625 :cvar INTERIOR_REDUCED: indicator that coefficient is defined on the
1193 caltinay 2169 interior of the PDE domain using a reduced
1194     integration order
1195 jfenwick 2625 :cvar BOUNDARY_REDUCED: indicator that coefficient is defined on the
1196 caltinay 2169 boundary of the PDE domain using a reduced
1197     integration order
1198 jfenwick 2625 :cvar CONTACT_REDUCED: indicator that coefficient is defined on the contact
1199 caltinay 2169 region within the PDE domain using a reduced
1200     integration order
1201 jfenwick 2625 :cvar SOLUTION: indicator that coefficient is defined through a solution of
1202 caltinay 2169 the PDE
1203 jfenwick 2625 :cvar REDUCED: indicator that coefficient is defined through a reduced
1204 caltinay 2169 solution of the PDE
1205 jfenwick 2625 :cvar BY_EQUATION: indicator that the dimension of the coefficient shape is
1206 caltinay 2169 defined by the number of PDE equations
1207 jfenwick 2625 :cvar BY_SOLUTION: indicator that the dimension of the coefficient shape is
1208 caltinay 2169 defined by the number of PDE solutions
1209 jfenwick 2625 :cvar BY_DIM: indicator that the dimension of the coefficient shape is
1210 caltinay 2169 defined by the spatial dimension
1211 jfenwick 2625 :cvar OPERATOR: indicator that the the coefficient alters the operator of
1212 caltinay 2169 the PDE
1213 jfenwick 2625 :cvar RIGHTHANDSIDE: indicator that the the coefficient alters the right
1214 caltinay 2169 hand side of the PDE
1215 jfenwick 2625 :cvar BOTH: indicator that the the coefficient alters the operator as well
1216 caltinay 2169 as the right hand side of the PDE
1217 jgs 149
1218 jgs 102 """
1219     INTERIOR=0
1220     BOUNDARY=1
1221     CONTACT=2
1222 jgs 149 SOLUTION=3
1223 jgs 150 REDUCED=4
1224 jgs 149 BY_EQUATION=5
1225     BY_SOLUTION=6
1226     BY_DIM=7
1227     OPERATOR=10
1228     RIGHTHANDSIDE=11
1229     BOTH=12
1230 gross 1072 INTERIOR_REDUCED=13
1231     BOUNDARY_REDUCED=14
1232     CONTACT_REDUCED=15
1233 jgs 149
1234 gross 1072 def __init__(self, where, pattern, altering):
1235 jgs 102 """
1236 caltinay 2169 Initialises a PDE coefficient type.
1237 jgs 151
1238 jfenwick 2625 :param where: describes where the coefficient lives
1239     :type where: one of `INTERIOR`, `BOUNDARY`, `CONTACT`, `SOLUTION`,
1240     `REDUCED`, `INTERIOR_REDUCED`, `BOUNDARY_REDUCED`,
1241     `CONTACT_REDUCED`
1242     :param pattern: describes the shape of the coefficient and how the shape
1243 caltinay 2169 is built for a given spatial dimension and numbers of
1244     equations and solutions in then PDE. For instance,
1245 jfenwick 2625 (`BY_EQUATION`,`BY_SOLUTION`,`BY_DIM`) describes a
1246 caltinay 2169 rank 3 coefficient which is instantiated as shape (3,2,2)
1247     in case of three equations and two solution components
1248     on a 2-dimensional domain. In the case of single equation
1249     and a single solution component the shape components
1250 jfenwick 2625 marked by `BY_EQUATION` or `BY_SOLUTION` are dropped.
1251 caltinay 2169 In this case the example would be read as (2,).
1252 jfenwick 2625 :type pattern: ``tuple`` of `BY_EQUATION`, `BY_SOLUTION`, `BY_DIM`
1253     :param altering: indicates what part of the PDE is altered if the
1254 caltinay 2169 coefficient is altered
1255 jfenwick 2625 :type altering: one of `OPERATOR`, `RIGHTHANDSIDE`, `BOTH`
1256 jgs 102 """
1257 gross 1841 super(PDECoef, self).__init__()
1258 jgs 102 self.what=where
1259     self.pattern=pattern
1260     self.altering=altering
1261 jgs 108 self.resetValue()
1262 jgs 102
1263 jgs 108 def resetValue(self):
1264     """
1265 caltinay 2169 Resets the coefficient value to the default.
1266 jgs 108 """
1267     self.value=escript.Data()
1268    
1269 jgs 149 def getFunctionSpace(self,domain,reducedEquationOrder=False,reducedSolutionOrder=False):
1270 jgs 102 """
1271 jfenwick 2625 Returns the `FunctionSpace` of the coefficient.
1272 jgs 102
1273 jfenwick 2625 :param domain: domain on which the PDE uses the coefficient
1274     :type domain: `Domain`
1275     :param reducedEquationOrder: True to indicate that reduced order is used
1276 caltinay 2169 to represent the equation
1277 jfenwick 2625 :type reducedEquationOrder: ``bool``
1278     :param reducedSolutionOrder: True to indicate that reduced order is used
1279 caltinay 2169 to represent the solution
1280 jfenwick 2625 :type reducedSolutionOrder: ``bool``
1281     :return: `FunctionSpace` of the coefficient
1282     :rtype: `FunctionSpace`
1283 jgs 102 """
1284 jgs 151 if self.what==self.INTERIOR:
1285 jgs 149 return escript.Function(domain)
1286 gross 1072 elif self.what==self.INTERIOR_REDUCED:
1287     return escript.ReducedFunction(domain)
1288 jgs 151 elif self.what==self.BOUNDARY:
1289 jgs 149 return escript.FunctionOnBoundary(domain)
1290 gross 1072 elif self.what==self.BOUNDARY_REDUCED:
1291     return escript.ReducedFunctionOnBoundary(domain)
1292 jgs 151 elif self.what==self.CONTACT:
1293 jgs 149 return escript.FunctionOnContactZero(domain)
1294 gross 1072 elif self.what==self.CONTACT_REDUCED:
1295     return escript.ReducedFunctionOnContactZero(domain)
1296 jgs 151 elif self.what==self.SOLUTION:
1297 jgs 149 if reducedEquationOrder and reducedSolutionOrder:
1298     return escript.ReducedSolution(domain)
1299     else:
1300     return escript.Solution(domain)
1301 jgs 151 elif self.what==self.REDUCED:
1302     return escript.ReducedSolution(domain)
1303 jgs 102
1304 jgs 108 def getValue(self):
1305     """
1306 caltinay 2169 Returns the value of the coefficient.
1307 jgs 149
1308 jfenwick 2625 :return: value of the coefficient
1309     :rtype: `Data`
1310 jgs 108 """
1311     return self.value
1312 jgs 148
1313 jgs 149 def setValue(self,domain,numEquations=1,numSolutions=1,reducedEquationOrder=False,reducedSolutionOrder=False,newValue=None):
1314 jgs 108 """
1315 caltinay 2169 Sets the value of the coefficient to a new value.
1316 jgs 149
1317 jfenwick 2625 :param domain: domain on which the PDE uses the coefficient
1318     :type domain: `Domain`
1319     :param numEquations: number of equations of the PDE
1320     :type numEquations: ``int``
1321     :param numSolutions: number of components of the PDE solution
1322     :type numSolutions: ``int``
1323     :param reducedEquationOrder: True to indicate that reduced order is used
1324 caltinay 2169 to represent the equation
1325 jfenwick 2625 :type reducedEquationOrder: ``bool``
1326     :param reducedSolutionOrder: True to indicate that reduced order is used
1327 caltinay 2169 to represent the solution
1328 jfenwick 2625 :type reducedSolutionOrder: ``bool``
1329     :param newValue: number of components of the PDE solution
1330     :type newValue: any object that can be converted into a
1331     `Data` object with the appropriate shape
1332     and `FunctionSpace`
1333     :raise IllegalCoefficientValue: if the shape of the assigned value does
1334 caltinay 2169 not match the shape of the coefficient
1335 jfenwick 2625 :raise IllegalCoefficientFunctionSpace: if unable to interpolate value
1336 caltinay 2169 to appropriate function space
1337 jgs 108 """
1338 jgs 148 if newValue==None:
1339     newValue=escript.Data()
1340     elif isinstance(newValue,escript.Data):
1341     if not newValue.isEmpty():
1342 gross 1072 if not newValue.getFunctionSpace() == self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder):
1343     try:
1344     newValue=escript.Data(newValue,self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder))
1345     except:
1346     raise IllegalCoefficientFunctionSpace,"Unable to interpolate coefficient to function space %s"%self.getFunctionSpace(domain)
1347 jgs 148 else:
1348 jgs 149 newValue=escript.Data(newValue,self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder))
1349 jgs 148 if not newValue.isEmpty():
1350     if not self.getShape(domain,numEquations,numSolutions)==newValue.getShape():
1351 jgs 149 raise IllegalCoefficientValue,"Expected shape of coefficient is %s but actual shape is %s."%(self.getShape(domain,numEquations,numSolutions),newValue.getShape())
1352 jgs 108 self.value=newValue
1353 jgs 148
1354 jgs 102 def isAlteringOperator(self):
1355     """
1356 caltinay 2169 Checks if the coefficient alters the operator of the PDE.
1357 jgs 149
1358 jfenwick 2625 :return: True if the operator of the PDE is changed when the
1359 caltinay 2169 coefficient is changed
1360 jfenwick 2625 :rtype: ``bool``
1361 caltinay 2169 """
1362 jgs 102 if self.altering==self.OPERATOR or self.altering==self.BOTH:
1363     return not None
1364     else:
1365     return None
1366    
1367     def isAlteringRightHandSide(self):
1368     """
1369 caltinay 2169 Checks if the coefficient alters the right hand side of the PDE.
1370 jgs 149
1371 jfenwick 2625 :rtype: ``bool``
1372     :return: True if the right hand side of the PDE is changed when the
1373     coefficient is changed, ``None`` otherwise.
1374 caltinay 2169 """
1375 jgs 102 if self.altering==self.RIGHTHANDSIDE or self.altering==self.BOTH:
1376     return not None
1377     else:
1378     return None
1379    
1380 jgs 148 def estimateNumEquationsAndNumSolutions(self,domain,shape=()):
1381 jgs 102 """
1382 caltinay 2169 Tries to estimate the number of equations and number of solutions if
1383     the coefficient has the given shape.
1384 jgs 102
1385 jfenwick 2625 :param domain: domain on which the PDE uses the coefficient
1386     :type domain: `Domain`
1387     :param shape: suggested shape of the coefficient
1388     :type shape: ``tuple`` of ``int`` values
1389     :return: the number of equations and number of solutions of the PDE if
1390 caltinay 2169 the coefficient has given shape. If no appropriate numbers
1391 jfenwick 2625 could be identified, ``None`` is returned
1392     :rtype: ``tuple`` of two ``int`` values or ``None``
1393 jgs 102 """
1394 jgs 148 dim=domain.getDim()
1395 jgs 102 if len(shape)>0:
1396     num=max(shape)+1
1397     else:
1398     num=1
1399     search=[]
1400 jgs 149 if self.definesNumEquation() and self.definesNumSolutions():
1401     for u in range(num):
1402     for e in range(num):
1403     search.append((e,u))
1404     search.sort(self.__CompTuple2)
1405     for item in search:
1406 jgs 148 s=self.getShape(domain,item[0],item[1])
1407 jgs 102 if len(s)==0 and len(shape)==0:
1408     return (1,1)
1409     else:
1410     if s==shape: return item
1411 jgs 149 elif self.definesNumEquation():
1412     for e in range(num,0,-1):
1413     s=self.getShape(domain,e,0)
1414     if len(s)==0 and len(shape)==0:
1415     return (1,None)
1416     else:
1417     if s==shape: return (e,None)
1418    
1419     elif self.definesNumSolutions():
1420     for u in range(num,0,-1):
1421     s=self.getShape(domain,0,u)
1422     if len(s)==0 and len(shape)==0:
1423     return (None,1)
1424     else:
1425     if s==shape: return (None,u)
1426 jgs 102 return None
1427 caltinay 2169
1428 jgs 149 def definesNumSolutions(self):
1429     """
1430 caltinay 2169 Checks if the coefficient allows to estimate the number of solution
1431     components.
1432 jgs 102
1433 jfenwick 2625 :return: True if the coefficient allows an estimate of the number of
1434 caltinay 2169 solution components, False otherwise
1435 jfenwick 2625 :rtype: ``bool``
1436 jgs 149 """
1437     for i in self.pattern:
1438     if i==self.BY_SOLUTION: return True
1439     return False
1440    
1441     def definesNumEquation(self):
1442     """
1443 caltinay 2169 Checks if the coefficient allows to estimate the number of equations.
1444 jgs 149
1445 jfenwick 2625 :return: True if the coefficient allows an estimate of the number of
1446 caltinay 2169 equations, False otherwise
1447 jfenwick 2625 :rtype: ``bool``
1448 jgs 149 """
1449     for i in self.pattern:
1450     if i==self.BY_EQUATION: return True
1451     return False
1452    
1453     def __CompTuple2(self,t1,t2):
1454     """
1455 caltinay 2169 Compares two tuples of possible number of equations and number of
1456     solutions.
1457 jgs 149
1458 jfenwick 2625 :param t1: the first tuple
1459     :param t2: the second tuple
1460     :return: 0, 1, or -1
1461 jgs 149 """
1462    
1463     dif=t1[0]+t1[1]-(t2[0]+t2[1])
1464     if dif<0: return 1
1465     elif dif>0: return -1
1466     else: return 0
1467    
1468 jgs 148 def getShape(self,domain,numEquations=1,numSolutions=1):
1469 jgs 149 """
1470 caltinay 2169 Builds the required shape of the coefficient.
1471 jgs 102
1472 jfenwick 2625 :param domain: domain on which the PDE uses the coefficient
1473     :type domain: `Domain`
1474     :param numEquations: number of equations of the PDE
1475     :type numEquations: ``int``
1476     :param numSolutions: number of components of the PDE solution
1477     :type numSolutions: ``int``
1478     :return: shape of the coefficient
1479     :rtype: ``tuple`` of ``int`` values
1480 jgs 149 """
1481     dim=domain.getDim()
1482     s=()
1483     for i in self.pattern:
1484     if i==self.BY_EQUATION:
1485 jgs 148 if numEquations>1: s=s+(numEquations,)
1486 jgs 149 elif i==self.BY_SOLUTION:
1487 jgs 148 if numSolutions>1: s=s+(numSolutions,)
1488 jgs 102 else:
1489     s=s+(dim,)
1490 jgs 149 return s
1491 jgs 102
1492 gross 2470 #====================================================================================================================
1493    
1494 gross 1841 class LinearProblem(object):
1495 jgs 102 """
1496 caltinay 2169 This is the base class to define a general linear PDE-type problem for
1497 jfenwick 2625 for an unknown function *u* on a given domain defined through a
1498     `Domain` object. The problem can be given as a single
1499 caltinay 2169 equation or as a system of equations.
1500 jgs 102
1501 caltinay 2169 The class assumes that some sort of assembling process is required to form
1502     a problem of the form
1503 jgs 151
1504 jfenwick 2625 *L u=f*
1505 jgs 102
1506 jfenwick 2625 where *L* is an operator and *f* is the right hand side. This operator
1507     problem will be solved to get the unknown *u*.
1508 caltinay 2169
1509 jgs 102 """
1510 jgs 148 def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):
1511 jgs 102 """
1512 caltinay 2169 Initializes a linear problem.
1513 jgs 102
1514 jfenwick 2625 :param domain: domain of the PDE
1515     :type domain: `Domain`
1516     :param numEquations: number of equations. If ``None`` the number of
1517 caltinay 2169 equations is extracted from the coefficients.
1518 jfenwick 2625 :param numSolutions: number of solution components. If ``None`` the number
1519 caltinay 2169 of solution components is extracted from the
1520     coefficients.
1521 jfenwick 2625 :param debug: if True debug information is printed
1522 jgs 148
1523 jgs 102 """
1524 gross 1841 super(LinearProblem, self).__init__()
1525 jgs 102
1526 jgs 148 self.__debug=debug
1527 jgs 104 self.__domain=domain
1528     self.__numEquations=numEquations
1529     self.__numSolutions=numSolutions
1530 gross 1841 self.__altered_coefficients=False
1531 jgs 149 self.__reduce_equation_order=False
1532     self.__reduce_solution_order=False
1533 gross 2474 self.__sym=False
1534     self.setSolverOptions()
1535 gross 1841 self.__is_RHS_valid=False
1536     self.__is_operator_valid=False
1537     self.__COEFFICIENTS={}
1538 gross 2474 self.__solution_rtol=1.e99
1539     self.__solution_atol=1.e99
1540     self.setSymmetryOff()
1541 gross 1841 # initialize things:
1542     self.resetAllCoefficients()
1543 gross 2474 self.initializeSystem()
1544 caltinay 2169 # ==========================================================================
1545 jgs 148 # general stuff:
1546 caltinay 2169 # ==========================================================================
1547 jgs 148 def __str__(self):
1548 jgs 149 """
1549 caltinay 2169 Returns a string representation of the PDE.
1550 jgs 149
1551 jfenwick 2625 :return: a simple representation of the PDE
1552     :rtype: ``str``
1553 jgs 149 """
1554 gross 1841 return "<LinearProblem %d>"%id(self)
1555 caltinay 2169 # ==========================================================================
1556 jgs 148 # debug :
1557 caltinay 2169 # ==========================================================================
1558 jgs 148 def setDebugOn(self):
1559 jgs 149 """
1560 caltinay 2169 Switches debug output on.
1561 jgs 108 """
1562 jgs 148 self.__debug=not None
1563    
1564     def setDebugOff(self):
1565 jgs 108 """
1566 caltinay 2169 Switches debug output off.
1567 jgs 148 """
1568     self.__debug=None
1569 jgs 108
1570 gross 2100 def setDebug(self, flag):
1571     """
1572 jfenwick 2625 Switches debug output on if ``flag`` is True otherwise it is switched off.
1573 caltinay 2169
1574 jfenwick 2625 :param flag: desired debug status
1575     :type flag: ``bool``
1576 gross 2100 """
1577     if flag:
1578     self.setDebugOn()
1579     else:
1580     self.setDebugOff()
1581    
1582 jgs 148 def trace(self,text):
1583     """
1584 caltinay 2169 Prints the text message if debug mode is switched on.
1585    
1586 jfenwick 2625 :param text: message to be printed
1587     :type text: ``string``
1588 jgs 102 """
1589 jgs 148 if self.__debug: print "%s: %s"%(str(self),text)
1590 jgs 102
1591 caltinay 2169 # ==========================================================================
1592 jgs 148 # some service functions:
1593 caltinay 2169 # ==========================================================================
1594 gross 1841 def introduceCoefficients(self,**coeff):
1595     """
1596 caltinay 2169 Introduces new coefficients into the problem.
1597 gross 1841
1598 jfenwick 2625 Use:
1599 gross 1841
1600 caltinay 2169 p.introduceCoefficients(A=PDECoef(...), B=PDECoef(...))
1601 gross 1841
1602 jfenwick 2625 to introduce the coefficients *A* and *B*.
1603 gross 1841 """
1604     for name, type in coeff.items():
1605     if not isinstance(type,PDECoef):
1606     raise ValueError,"coefficient %s has no type."%name
1607     self.__COEFFICIENTS[name]=type
1608     self.__COEFFICIENTS[name].resetValue()
1609     self.trace("coefficient %s has been introduced."%name)
1610    
1611 jgs 148 def getDomain(self):
1612 jgs 102 """
1613 caltinay 2169 Returns the domain of the PDE.
1614 jgs 102
1615 jfenwick 2625 :return: the domain of the PDE
1616     :rtype: `Domain`
1617 jgs 108 """
1618 jgs 148 return self.__domain
1619 gross 2535 def getDomainStatus(self):
1620     """
1621     Return the status indicator of the domain
1622     """
1623     return self.getDomain().getStatus()
1624 jgs 122
1625 gross 2535 def getSystemStatus(self):
1626     """
1627     Return the domain status used to build the current system
1628     """
1629     return self.__system_status
1630     def setSystemStatus(self,status=None):
1631     """
1632 jfenwick 2625 Sets the system status to ``status`` if ``status`` is not present the
1633 gross 2535 current status of the domain is used.
1634     """
1635     if status == None:
1636     self.__system_status=self.getDomainStatus()
1637     else:
1638     self.__system_status=status
1639    
1640 jgs 148 def getDim(self):
1641 jgs 108 """
1642 caltinay 2169 Returns the spatial dimension of the PDE.
1643 jgs 108
1644 jfenwick 2625 :return: the spatial dimension of the PDE domain
1645     :rtype: ``int``
1646 jgs 148 """
1647     return self.getDomain().getDim()
1648 jgs 102
1649 jgs 148 def getNumEquations(self):
1650     """
1651 caltinay 2169 Returns the number of equations.
1652 jgs 102
1653 jfenwick 2625 :return: the number of equations
1654     :rtype: ``int``
1655     :raise UndefinedPDEError: if the number of equations is not specified yet
1656 jgs 148 """
1657     if self.__numEquations==None:
1658 gross 1859 if self.__numSolutions==None:
1659     raise UndefinedPDEError,"Number of equations is undefined. Please specify argument numEquations."
1660     else:
1661     self.__numEquations=self.__numSolutions
1662     return self.__numEquations
1663 jgs 147
1664 jgs 148 def getNumSolutions(self):
1665     """
1666 caltinay 2169 Returns the number of unknowns.
1667 jgs 147
1668 jfenwick 2625 :return: the number of unknowns
1669     :rtype: ``int``
1670     :raise UndefinedPDEError: if the number of unknowns is not specified yet
1671 jgs 148 """
1672     if self.__numSolutions==None:
1673 gross 1859 if self.__numEquations==None:
1674     raise UndefinedPDEError,"Number of solution is undefined. Please specify argument numSolutions."
1675     else:
1676     self.__numSolutions=self.__numEquations
1677     return self.__numSolutions
1678 jgs 148
1679 jgs 149 def reduceEquationOrder(self):
1680     """
1681 caltinay 2169 Returns the status of order reduction for the equation.
1682 jgs 149
1683 jfenwick 2625 :return: True if reduced interpolation order is used for the
1684 caltinay 2169 representation of the equation, False otherwise
1685 jfenwick 2625 :rtype: `bool`
1686 jgs 149 """
1687     return self.__reduce_equation_order
1688    
1689     def reduceSolutionOrder(self):
1690     """
1691 caltinay 2169 Returns the status of order reduction for the solution.
1692 jgs 149
1693 jfenwick 2625 :return: True if reduced interpolation order is used for the
1694 caltinay 2169 representation of the solution, False otherwise
1695 jfenwick 2625 :rtype: `bool`
1696 jgs 149 """
1697     return self.__reduce_solution_order
1698 jgs 151
1699 jgs 108 def getFunctionSpaceForEquation(self):
1700     """
1701 jfenwick 2625 Returns the `FunctionSpace` used to discretize
1702 caltinay 2169 the equation.
1703 jgs 148
1704 jfenwick 2625 :return: representation space of equation
1705     :rtype: `FunctionSpace`
1706 jgs 108 """
1707 jgs 149 if self.reduceEquationOrder():
1708     return escript.ReducedSolution(self.getDomain())
1709     else:
1710     return escript.Solution(self.getDomain())
1711 jgs 108
1712     def getFunctionSpaceForSolution(self):
1713     """
1714 jfenwick 2625 Returns the `FunctionSpace` used to represent
1715 caltinay 2169 the solution.
1716 jgs 148
1717 jfenwick 2625 :return: representation space of solution
1718     :rtype: `FunctionSpace`
1719 jgs 108 """
1720 jgs 149 if self.reduceSolutionOrder():
1721     return escript.ReducedSolution(self.getDomain())
1722     else:
1723     return escript.Solution(self.getDomain())
1724 jgs 108
1725 caltinay 2169 # ==========================================================================
1726 jgs 148 # solver settings:
1727 caltinay 2169 # ==========================================================================
1728 gross 2474 def setSolverOptions(self,options=None):
1729 jgs 102 """
1730 gross 2474 Sets the solver options.
1731 jgs 148
1732 jfenwick 2625 :param options: the new solver options. If equal ``None``, the solver options are set to the default.
1733     :type options: `SolverOptions` or ``None``
1734     :note: The symmetry flag of options is overwritten by the symmetry flag of the `LinearProblem`.
1735 jgs 102 """
1736 gross 2474 if options==None:
1737     self.__solver_options=SolverOptions()
1738     elif isinstance(options, SolverOptions):
1739     self.__solver_options=options
1740     else:
1741     raise ValueError,"options must be a SolverOptions object."
1742     self.__solver_options.setSymmetry(self.__sym)
1743    
1744     def getSolverOptions(self):
1745 jgs 148 """
1746 gross 2474 Returns the solver options
1747    
1748 jfenwick 2625 :rtype: `SolverOptions`
1749 jgs 148 """
1750 gross 2474 self.__solver_options.setSymmetry(self.__sym)
1751     return self.__solver_options
1752    
1753 jgs 148 def isUsingLumping(self):
1754     """
1755 caltinay 2169 Checks if matrix lumping is the current solver method.
1756 jgs 148
1757 jfenwick 2625 :return: True if the current solver method is lumping
1758     :rtype: ``bool``
1759 jgs 148 """
1760 gross 3379 return self.getSolverOptions().getSolverMethod() in [ SolverOptions.ROWSUM_LUMPING, SolverOptions.HRZ_LUMPING ]
1761 caltinay 2169 # ==========================================================================
1762 jgs 148 # symmetry flag:
1763 caltinay 2169 # ==========================================================================
1764     def isSymmetric(self):
1765 jgs 102 """
1766 caltinay 2169 Checks if symmetry is indicated.
1767 jgs 149
1768 jfenwick 2625 :return: True if a symmetric PDE is indicated, False otherwise
1769     :rtype: ``bool``
1770     :note: the method is equivalent to use getSolverOptions().isSymmetric()
1771 jgs 102 """
1772 gross 2474 self.getSolverOptions().isSymmetric()
1773 jgs 102
1774 caltinay 2169 def setSymmetryOn(self):
1775 jgs 102 """
1776 gross 2474 Sets the symmetry flag.
1777 jfenwick 2625 :note: The method overwrites the symmetry flag set by the solver options
1778 jgs 102 """
1779 gross 2474 self.__sym=True
1780     self.getSolverOptions().setSymmetryOn()
1781 jgs 102
1782 caltinay 2169 def setSymmetryOff(self):
1783 jgs 102 """
1784 caltinay 2169 Clears the symmetry flag.
1785 jfenwick 2625 :note: The method overwrites the symmetry flag set by the solver options
1786 jgs 102 """
1787 gross 2474 self.__sym=False
1788     self.getSolverOptions().setSymmetryOff()
1789 jgs 102
1790 gross 2474 def setSymmetry(self,flag=False):
1791 jgs 148 """
1792 jfenwick 2625 Sets the symmetry flag to ``flag``.
1793 jgs 149
1794 jfenwick 2625 :param flag: If True, the symmetry flag is set otherwise reset.
1795     :type flag: ``bool``
1796     :note: The method overwrites the symmetry flag set by the solver options
1797 jgs 148 """
1798 gross 2474 self.getSolverOptions().setSymmetry(flag)
1799 caltinay 2169 # ==========================================================================
1800 jgs 148 # function space handling for the equation as well as the solution
1801 caltinay 2169 # ==========================================================================
1802     def setReducedOrderOn(self):
1803 jgs 102 """
1804 caltinay 2169 Switches reduced order on for solution and equation representation.
1805 jgs 149
1806 jfenwick 2625 :raise RuntimeError: if order reduction is altered after a coefficient has
1807 caltinay 2169 been set
1808 jgs 102 """
1809     self.setReducedOrderForSolutionOn()
1810     self.setReducedOrderForEquationOn()
1811    
1812 caltinay 2169 def setReducedOrderOff(self):
1813 jgs 102 """
1814 caltinay 2169 Switches reduced order off for solution and equation representation
1815 jgs 149
1816 jfenwick 2625 :raise RuntimeError: if order reduction is altered after a coefficient has
1817 caltinay 2169 been set
1818 jgs 102 """
1819     self.setReducedOrderForSolutionOff()
1820     self.setReducedOrderForEquationOff()
1821    
1822 caltinay 2169 def setReducedOrderTo(self,flag=False):
1823 jgs 102 """
1824 caltinay 2169 Sets order reduction state for both solution and equation representation
1825     according to flag.
1826    
1827 jfenwick 2625 :param flag: if True, the order reduction is switched on for both solution
1828 caltinay 2169 and equation representation, otherwise or if flag is not
1829     present order reduction is switched off
1830 jfenwick 2625 :type flag: ``bool``
1831     :raise RuntimeError: if order reduction is altered after a coefficient has
1832 caltinay 2169 been set
1833 jgs 102 """
1834     self.setReducedOrderForSolutionTo(flag)
1835     self.setReducedOrderForEquationTo(flag)
1836    
1837 jgs 148
1838 caltinay 2169 def setReducedOrderForSolutionOn(self):
1839 jgs 102 """
1840 caltinay 2169 Switches reduced order on for solution representation.
1841 jgs 149
1842 jfenwick 2625 :raise RuntimeError: if order reduction is altered after a coefficient has
1843 caltinay 2169 been set
1844 jgs 102 """
1845 jgs 149 if not self.__reduce_solution_order:
1846     if self.__altered_coefficients:
1847     raise RuntimeError,"order cannot be altered after coefficients have been defined."
1848 caltinay 2169 self.trace("Reduced order is used for solution representation.")
1849 jgs 149 self.__reduce_solution_order=True
1850 gross 1841 self.initializeSystem()
1851 jgs 102
1852 caltinay 2169 def setReducedOrderForSolutionOff(self):
1853 jgs 102 """
1854 caltinay 2169 Switches reduced order off for solution representation
1855 jgs 149
1856 jfenwick 2625 :raise RuntimeError: if order reduction is altered after a coefficient has
1857 caltinay 2169 been set.
1858 jgs 102 """
1859 jgs 149 if self.__reduce_solution_order:
1860     if self.__altered_coefficients:
1861     raise RuntimeError,"order cannot be altered after coefficients have been defined."
1862 jgs 148 self.trace("Full order is used to interpolate solution.")
1863 jgs 149 self.__reduce_solution_order=False
1864 gross 1841 self.initializeSystem()
1865 jgs 102
1866 caltinay 2169 def setReducedOrderForSolutionTo(self,flag=False):
1867 jgs 102 """
1868 caltinay 2169 Sets order reduction state for solution representation according to flag.
1869 jgs 102
1870 jfenwick 2625 :param flag: if flag is True, the order reduction is switched on for
1871 caltinay 2169 solution representation, otherwise or if flag is not present
1872     order reduction is switched off
1873 jfenwick 2625 :type flag: ``bool``
1874     :raise RuntimeError: if order reduction is altered after a coefficient has
1875 caltinay 2169 been set
1876 jgs 102 """
1877     if flag:
1878     self.setReducedOrderForSolutionOn()
1879     else:
1880     self.setReducedOrderForSolutionOff()
1881 jgs 148
1882 caltinay 2169 def setReducedOrderForEquationOn(self):
1883 jgs 102 """
1884 caltinay 2169 Switches reduced order on for equation representation.
1885 jgs 149
1886 jfenwick 2625 :raise RuntimeError: if order reduction is altered after a coefficient has
1887 caltinay 2169 been set
1888 jgs 102 """
1889 jgs 149 if not self.__reduce_equation_order:
1890     if self.__altered_coefficients:
1891     raise RuntimeError,"order cannot be altered after coefficients have been defined."
1892 jgs 148 self.trace("Reduced order is used for test functions.")
1893 jgs 149 self.__reduce_equation_order=True
1894 gross 1841 self.initializeSystem()
1895 jgs 102
1896 caltinay 2169 def setReducedOrderForEquationOff(self):
1897 jgs 102 """
1898 caltinay 2169 Switches reduced order off for equation representation.
1899 jgs 149
1900 jfenwick 2625 :raise RuntimeError: if order reduction is altered after a coefficient has
1901 caltinay 2169 been set
1902 jgs 102 """
1903 jgs 149 if self.__reduce_equation_order:
1904     if self.__altered_coefficients:
1905     raise RuntimeError,"order cannot be altered after coefficients have been defined."
1906 jgs 148 self.trace("Full order is used for test functions.")
1907 jgs 149 self.__reduce_equation_order=False
1908 gross 1841 self.initializeSystem()
1909 jgs 102
1910 caltinay 2169 def setReducedOrderForEquationTo(self,flag=False):
1911 jgs 102 """
1912 caltinay 2169 Sets order reduction state for equation representation according to flag.
1913 jgs 102
1914 jfenwick 2625 :param flag: if flag is True, the order reduction is switched on for
1915 caltinay 2169 equation representation, otherwise or if flag is not present
1916     order reduction is switched off
1917 jfenwick 2625 :type flag: ``bool``
1918     :raise RuntimeError: if order reduction is altered after a coefficient has
1919 caltinay 2169 been set
1920 jgs 102 """
1921     if flag:
1922     self.setReducedOrderForEquationOn()
1923     else:
1924     self.setReducedOrderForEquationOff()
1925 gross 2474 def getOperatorType(self):
1926 gross 1841 """
1927 caltinay 2169 Returns the current system type.
1928 gross 1841 """
1929 gross 2474 return self.__operator_type
1930 jgs 148
1931 gross 1841 def checkSymmetricTensor(self,name,verbose=True):
1932     """
1933 caltinay 2169 Tests a coefficient for symmetry.
1934 jgs 148
1935 jfenwick 2625 :param name: name of the coefficient
1936     :type name: ``str``
1937     :param verbose: if set to True or not present a report on coefficients
1938 caltinay 2169 which break the symmetry is printed.
1939 jfenwick 2625 :type verbose: ``bool``
1940     :return: True if coefficient ``name`` is symmetric
1941     :rtype: ``bool``
1942 gross 1841 """
1943 gross 2474 SMALL_TOLERANCE=util.EPSILON*10.
1944 gross 1841 A=self.getCoefficient(name)
1945     verbose=verbose or self.__debug
1946     out=True
1947     if not A.isEmpty():
1948 gross 2474 tol=util.Lsup(A)*SMALL_TOLERANCE
1949 gross 1841 s=A.getShape()
1950     if A.getRank() == 4:
1951     if s[0]==s[2] and s[1] == s[3]:
1952     for i in range(s[0]):
1953     for j in range(s[1]):
1954     for k in range(s[2]):
1955     for l in range(s[3]):
1956     if util.Lsup(A[i,j,k,l]-A[k,l,i,j])>tol:
1957     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)
1958     out=False
1959     else:
1960     if verbose: print "non-symmetric problem because of inappropriate shape %s of coefficient %s."%(s,name)
1961     out=False
1962     elif A.getRank() == 2:
1963     if s[0]==s[1]:
1964     for j in range(s[0]):
1965     for l in range(s[1]):
1966     if util.Lsup(A[j,l]-A[l,j])>tol:
1967     if verbose: print "non-symmetric problem because %s[%d,%d]!=%s[%d,%d]"%(name,j,l,name,l,j)
1968     out=False
1969     else:
1970     if verbose: print "non-symmetric problem because of inappropriate shape %s of coefficient %s."%(s,name)
1971     out=False
1972     elif A.getRank() == 0:
1973     pass
1974     else:
1975     raise ValueError,"Cannot check rank %s of %s."%(A.getRank(),name)
1976     return out
1977 jgs 148
1978 gross 1841 def checkReciprocalSymmetry(self,name0,name1,verbose=True):
1979     """
1980 caltinay 2169 Tests two coefficients for reciprocal symmetry.
1981 jgs 148
1982 jfenwick 2625 :param name0: name of the first coefficient
1983     :type name0: ``str``
1984     :param name1: name of the second coefficient
1985     :type name1: ``str``
1986     :param verbose: if set to True or not present a report on coefficients
1987 caltinay 2169 which break the symmetry is printed
1988 jfenwick 2625 :type verbose: ``bool``
1989     :return: True if coefficients ``name0`` and ``name1`` are reciprocally
1990 caltinay 2169 symmetric.
1991 jfenwick 2625 :rtype: ``bool``
1992 gross 1841 """
1993 gross 2474 SMALL_TOLERANCE=util.EPSILON*10.
1994 gross 1841 B=self.getCoefficient(name0)
1995     C=self.getCoefficient(name1)
1996     verbose=verbose or self.__debug
1997     out=True
1998     if B.isEmpty() and not C.isEmpty():
1999     if verbose: print "non-symmetric problem because %s is not present but %s is"%(name0,name1)
2000     out=False
2001     elif not B.isEmpty() and C.isEmpty():
2002     if verbose: print "non-symmetric problem because %s is not present but %s is"%(name0,name1)
2003     out=False
2004     elif not B.isEmpty() and not C.isEmpty():
2005     sB=B.getShape()
2006     sC=C.getShape()
2007 gross 2474 tol=(util.Lsup(B)+util.Lsup(C))*SMALL_TOLERANCE/2.
2008 gross 1841 if len(sB) != len(sC):
2009     if verbose: print "non-symmetric problem because ranks of %s (=%s) and %s (=%s) are different."%(name0,len(sB),name1,len(sC))
2010     out=False
2011     else:
2012     if len(sB)==0:
2013     if util.Lsup(B-C)>tol:
2014     if verbose: print "non-symmetric problem because %s!=%s"%(name0,name1)
2015     out=False
2016     elif len(sB)==1:
2017     if sB[0]==sC[0]:
2018     for j in range(sB[0]):
2019     if util.Lsup(B[j]-C[j])>tol:
2020     if verbose: print "non-symmetric PDE because %s[%d]!=%s[%d]"%(name0,j,name1,j)
2021     out=False
2022     else:
2023     if verbose: print "non-symmetric problem because of inappropriate shapes %s and %s of coefficients %s and %s, respectively."%(sB,sC,name0,name1)
2024     elif len(sB)==3:
2025     if sB[0]==sC[1] and sB[1]==sC[2] and sB[2]==sC[0]:
2026     for i in range(sB[0]):
2027     for j in range(sB[1]):
2028     for k in range(sB[2]):
2029     if util.Lsup(B[i,j,k]-C[k,i,j])>tol:
2030 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)
2031 gross 1841 out=False
2032     else:
2033     if verbose: print "non-symmetric problem because of inappropriate shapes %s and %s of coefficients %s and %s, respectively."%(sB,sC,name0,name1)
2034 jgs 148 else:
2035 gross 1841 raise ValueError,"Cannot check rank %s of %s and %s."%(len(sB),name0,name1)
2036     return out
2037 jgs 148
2038 caltinay 2169 def getCoefficient(self,name):
2039 jgs 108 """
2040 jfenwick 2625 Returns the value of the coefficient ``name``.
2041 jgs 108
2042 jfenwick 2625 :param name: name of the coefficient requested
2043     :type name: ``string``
2044     :return: the value of the coefficient
2045     :rtype: `Data`
2046     :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE
2047 jgs 148 """
2048     if self.hasCoefficient(name):
2049 gross 1841 return self.__COEFFICIENTS[name].getValue()
2050 jgs 148 else:
2051     raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
2052 jgs 108
2053 caltinay 2169 def hasCoefficient(self,name):
2054 jgs 148 """
2055 jfenwick 2625 Returns True if ``name`` is the name of a coefficient.
2056 jgs 108
2057 jfenwick 2625 :param name: name of the coefficient enquired
2058     :type name: ``string``
2059     :return: True if ``name`` is the name of a coefficient of the general PDE,
2060 caltinay 2169 False otherwise
2061 jfenwick 2625 :rtype: ``bool``
2062 jgs 148 """
2063 gross 1841 return self.__COEFFICIENTS.has_key(name)
2064 jgs 108
2065 caltinay 2169 def createCoefficient(self, name):
2066 jgs 148 """
2067 jfenwick 2625 Creates a `Data` object corresponding to coefficient
2068     ``name``.
2069 jgs 108
2070 jfenwick 2625 :return: the coefficient ``name`` initialized to 0
2071     :rtype: `Data`
2072     :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE
2073 jgs 148 """
2074     if self.hasCoefficient(name):
2075     return escript.Data(0.,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name))
2076     else:
2077     raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
2078    
2079 caltinay 2169 def getFunctionSpaceForCoefficient(self,name):
2080 jgs 148 """
2081 jfenwick 2625 Returns the `FunctionSpace` to be used for
2082     coefficient ``name``.
2083 jgs 148
2084 jfenwick 2625 :param name: name of the coefficient enquired
2085     :type name: ``string``
2086     :return: the function space to be used for coefficient ``name``
2087     :rtype: `FunctionSpace`
2088     :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE
2089 jgs 148 """
2090     if self.hasCoefficient(name):
2091 gross 1841 return self.__COEFFICIENTS[name].getFunctionSpace(self.getDomain())
2092 jgs 148 else:
2093     raise ValueError,"unknown coefficient %s requested"%name
2094 gross 1841
2095 caltinay 2169 def getShapeOfCoefficient(self,name):
2096 jgs 148 """
2097 jfenwick 2625 Returns the shape of the coefficient ``name``.
2098 jgs 148
2099 jfenwick 2625 :param name: name of the coefficient enquired
2100     :type name: ``string``
2101     :return: the shape of the coefficient ``name``
2102     :rtype: ``tuple`` of ``int``
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].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())
2107 jgs 148 else:
2108     raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
2109    
2110 caltinay 2169 def resetAllCoefficients(self):
2111 jgs 148 """
2112 caltinay 2169 Resets all coefficients to their default values.
2113 jgs 148 """
2114 gross 1841 for i in self.__COEFFICIENTS.iterkeys():
2115     self.__COEFFICIENTS[i].resetValue()
2116 jgs 148
2117 caltinay 2169 def alteredCoefficient(self,name):
2118 jgs 148 """
2119 jfenwick 2625 Announces that coefficient ``name`` has been changed.
2120 jgs 148
2121 jfenwick 2625 :param name: name of the coefficient affected
2122     :type name: ``string``
2123     :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE
2124     :note: if ``name`` is q or r, the method will not trigger a rebuild of the
2125 caltinay 2169 system as constraints are applied to the solved system.
2126 jgs 148 """
2127     if self.hasCoefficient(name):
2128     self.trace("Coefficient %s has been altered."%name)
2129 jgs 149 if not ((name=="q" or name=="r") and self.isUsingLumping()):
2130 gross 1841 if self.__COEFFICIENTS[name].isAlteringOperator(): self.invalidateOperator()
2131     if self.__COEFFICIENTS[name].isAlteringRightHandSide(): self.invalidateRightHandSide()
2132 jgs 148 else:
2133     raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
2134    
2135 caltinay 2169 def validSolution(self):
2136 gross 1841 """
2137 caltinay 2169 Marks the solution as valid.
2138 gross 1841 """
2139     self.__is_solution_valid=True
2140 jgs 149
2141 caltinay 2169 def invalidateSolution(self):
2142 gross 1841 """
2143 caltinay 2169 Indicates the PDE has to be resolved if the solution is requested.
2144 gross 1841 """
2145     self.trace("System will be resolved.")
2146     self.__is_solution_valid=False
2147 jgs 148
2148 gross 1841 def isSolutionValid(self):
2149     """
2150 caltinay 2169 Returns True if the solution is still valid.
2151 gross 1841 """
2152 gross 2535 if not self.getDomainStatus()==self.getSystemStatus(): self.invalidateSolution()
2153 gross 2474 if self.__solution_rtol>self.getSolverOptions().getTolerance() or \
2154     self.__solution_atol>self.getSolverOptions().getAbsoluteTolerance():
2155     self.invalidateSolution()
2156 gross 1841 return self.__is_solution_valid
2157    
2158     def validOperator(self):
2159     """
2160 caltinay 2169 Marks the operator as valid.
2161 gross 1841 """
2162     self.__is_operator_valid=True
2163    
2164 caltinay 2169 def invalidateOperator(self):
2165 gross 1841 """
2166 caltinay 2169 Indicates the operator has to be rebuilt next time it is used.
2167 gross 1841 """
2168     self.trace("Operator will be rebuilt.")
2169     self.invalidateSolution()
2170     self.__is_operator_valid=False
2171    
2172     def isOperatorValid(self):
2173     """
2174 caltinay 2169 Returns True if the operator is still valid.
2175 gross 1841 """
2176 gross 2535 if not self.getDomainStatus()==self.getSystemStatus(): self.invalidateOperator()
2177     if not self.getRequiredOperatorType()==self.getOperatorType(): self.invalidateOperator()
2178 gross 1841 return self.__is_operator_valid
2179    
2180     def validRightHandSide(self):
2181     """
2182 caltinay 2169 Marks the right hand side as valid.
2183 gross 1841 """
2184     self.__is_RHS_valid=True
2185    
2186 caltinay 2169 def invalidateRightHandSide(self):
2187 gross 1841 """
2188 caltinay 2169 Indicates the right hand side has to be rebuilt next time it is used.
2189 gross 1841 """
2190 gross 2535 self.trace("Right hand side has to be rebuilt.")
2191 gross 1841 self.invalidateSolution()
2192     self.__is_RHS_valid=False
2193    
2194     def isRightHandSideValid(self):
2195     """
2196 caltinay 2169 Returns True if the operator is still valid.
2197 gross 1841 """
2198 gross 2535 if not self.getDomainStatus()==self.getSystemStatus(): self.invalidateRightHandSide()
2199 gross 1841 return self.__is_RHS_valid
2200    
2201 caltinay 2169 def invalidateSystem(self):
2202 gross 1841 """
2203 caltinay 2169 Announces that everything has to be rebuilt.
2204 gross 1841 """
2205     self.invalidateSolution()
2206     self.invalidateOperator()
2207     self.invalidateRightHandSide()
2208    
2209     def isSystemValid(self):
2210     """
2211 caltinay 2169 Returns True if the system (including solution) is still vaild.
2212 gross 1841 """
2213     return self.isSolutionValid() and self.isOperatorValid() and self.isRightHandSideValid()
2214    
2215 caltinay 2169 def initializeSystem(self):
2216 gross 1841 """
2217 caltinay 2169 Resets the system clearing the operator, right hand side and solution.
2218 gross 1841 """
2219     self.trace("New System has been created.")
2220 gross 2474 self.__operator_type=None
2221 gross 2535 self.setSystemStatus()
2222 gross 1841 self.__operator=escript.Operator()
2223     self.__righthandside=escript.Data()
2224     self.__solution=escript.Data()
2225     self.invalidateSystem()
2226    
2227 caltinay 2169 def getOperator(self):
2228 gross 1841 """
2229 caltinay 2169 Returns the operator of the linear problem.
2230 gross 1841
2231 jfenwick 2625 :return: the operator of the problem
2232 gross 1841 """
2233     return self.getSystem()[0]
2234    
2235 caltinay 2169 def getRightHandSide(self):
2236 gross 1841 """
2237 caltinay 2169 Returns the right hand side of the linear problem.
2238 gross 1841
2239 jfenwick 2625 :return: the right hand side of the problem
2240     :rtype: `Data`
2241 gross 1841 """
2242     return self.getSystem()[1]
2243    
2244 caltinay 2169 def createRightHandSide(self):
2245 gross 1841 """
2246 caltinay 2169 Returns an instance of a new right hand side.
2247 gross 1841 """
2248     self.trace("New right hand side is allocated.")
2249     if self.getNumEquations()>1:
2250     return escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True)
2251     else:
2252     return escript.Data(0.,(),self.getFunctionSpaceForEquation(),True)
2253    
2254 caltinay 2169 def createSolution(self):
2255 gross 1841 """
2256 caltinay 2169 Returns an instance of a new solution.
2257 gross 1841 """
2258     self.trace("New solution is allocated.")
2259     if self.getNumSolutions()>1:
2260     return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)
2261     else:
2262     return escript.Data(0.,(),self.getFunctionSpaceForSolution(),True)
2263    
2264 caltinay 2169 def resetSolution(self):
2265 gross 1841 """
2266 caltinay 2169 Sets the solution to zero.
2267 gross 1841 """
2268     if self.__solution.isEmpty():
2269     self.__solution=self.createSolution()
2270     else:
2271     self.__solution.setToZero()
2272     self.trace("Solution is reset to zero.")
2273 caltinay 2169
2274 gross 2987 def setSolution(self,u, validate=True):
2275 gross 1841 """
2276 gross 2474 Sets the solution assuming that makes the system valid with the tolrance
2277     defined by the solver options
2278 gross 1841 """
2279 gross 2987 if validate:
2280     self.__solution_rtol=self.getSolverOptions().getTolerance()
2281     self.__solution_atol=self.getSolverOptions().getAbsoluteTolerance()
2282     self.validSolution()
2283 gross 1841 self.__solution=u
2284     def getCurrentSolution(self):
2285     """
2286 caltinay 2169 Returns the solution in its current state.
2287 gross 1841 """
2288     if self.__solution.isEmpty(): self.__solution=self.createSolution()
2289     return self.__solution
2290    
2291 caltinay 2169 def resetRightHandSide(self):
2292 gross 1841 """
2293 caltinay 2169 Sets the right hand side to zero.
2294 gross 1841 """
2295     if self.__righthandside.isEmpty():
2296     self.__righthandside=self.createRightHandSide()
2297     else:
2298     self.__righthandside.setToZero()
2299     self.trace("Right hand side is reset to zero.")
2300    
2301     def getCurrentRightHandSide(self):
2302     """
2303 caltinay 2169 Returns the right hand side in its current state.
2304 gross 1841 """
2305     if self.__righthandside.isEmpty(): self.__righthandside=self.createRightHandSide()
2306     return self.__righthandside
2307    
2308 caltinay 2169 def resetOperator(self):
2309 gross 1841 """
2310 caltinay 2169 Makes sure that the operator is instantiated and returns it initialized
2311     with zeros.
2312 gross 1841 """
2313 gross 2474 if self.getOperatorType() == None:
2314 gross 1841 if self.isUsingLumping():
2315     self.__operator=self.createSolution()
2316     else:
2317     self.__operator=self.createOperator()
2318 gross 2474 self.__operator_type=self.getRequiredOperatorType()
2319 gross 1841 else:
2320     if self.isUsingLumping():
2321     self.__operator.setToZero()
2322     else:
2323 gross 2551 if self.getOperatorType() == self.getRequiredOperatorType():
2324     self.__operator.resetValues()
2325     else:
2326     self.__operator=self.createOperator()
2327     self.__operator_type=self.getRequiredOperatorType()
2328 gross 1841 self.trace("Operator reset to zero")
2329    
2330     def getCurrentOperator(self):
2331     """
2332 caltinay 2169 Returns the operator in its current state.
2333 gross 1841 """
2334     return self.__operator
2335    
2336 caltinay 2169 def setValue(self,**coefficients):
2337 jgs 148 """
2338 caltinay 2169 Sets new values to coefficients.
2339 jgs 148
2340 jfenwick 2625 :raise IllegalCoefficient: if an unknown coefficient keyword is used
2341 jgs 148 """
2342 jgs 108 # check if the coefficients are legal:
2343     for i in coefficients.iterkeys():
2344     if not self.hasCoefficient(i):
2345 jgs 148 raise IllegalCoefficient,"Attempt to set unknown coefficient %s"%i
2346 jgs 108 # if the number of unknowns or equations is still unknown we try to estimate them:
2347 jgs 148 if self.__numEquations==None or self.__numSolutions==None:
2348 jgs 108 for i,d in coefficients.iteritems():
2349     if hasattr(d,"shape"):
2350     s=d.shape
2351     elif hasattr(d,"getShape"):
2352     s=d.getShape()
2353     else:
2354 jfenwick 2455 s=numpy.array(d).shape
2355 jgs 108 if s!=None:
2356     # get number of equations and number of unknowns:
2357 gross 1841 res=self.__COEFFICIENTS[i].estimateNumEquationsAndNumSolutions(self.getDomain(),s)
2358 jgs 108 if res==None:
2359 jgs 148 raise IllegalCoefficientValue,"Illegal shape %s of coefficient %s"%(s,i)
2360 jgs 108 else:
2361 jgs 148 if self.__numEquations==None: self.__numEquations=res[0]
2362     if self.__numSolutions==None: self.__numSolutions=res[1]
2363 caltinay 2169 if self.__numEquations==None: raise UndefinedPDEError,"unidentified number of equations"
2364     if self.__numSolutions==None: raise UndefinedPDEError,"unidentified number of solutions"
2365 jgs 108 # now we check the shape of the coefficient if numEquations and numSolutions are set:
2366     for i,d in coefficients.iteritems():
2367 jgs 148 try:
2368 gross 1841 self.__COEFFICIENTS[i].setValue(self.getDomain(),
2369 caltinay 2169 self.getNumEquations(),self.getNumSolutions(),
2370     self.reduceEquationOrder(),self.reduceSolutionOrder(),d)
2371 gross 1072 self.alteredCoefficient(i)
2372     except IllegalCoefficientFunctionSpace,m:
2373     # if the function space is wrong then we try the reduced version:
2374     i_red=i+"_reduced"
2375 gross 1841 if (not i_red in coefficients.keys()) and i_red in self.__COEFFICIENTS.keys():
2376 gross 1072 try:
2377 gross 1841 self.__COEFFICIENTS[i_red].setValue(self.getDomain(),
2378 gross