/[escript]/trunk/escript/py_src/linearPDEs.py
ViewVC logotype

Annotation of /trunk/escript/py_src/linearPDEs.py

Parent Directory Parent Directory | Revision Log Revision Log


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