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