/[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 3471 - (hide annotations)
Tue Mar 15 04:23:54 2011 UTC (8 years, 8 months ago) by caltinay
File MIME type: text/x-python
File size: 161133 byte(s)
Branch to investigate how to incorporate sympy or other symbolic evaluation
libraries into escript.

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