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