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