/[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 3094 - (hide annotations)
Fri Aug 13 08:38:06 2010 UTC (9 years, 1 month ago) by gross
File MIME type: text/x-python
File size: 153310 byte(s)
The MPI and sequational GAUSS_SEIDEL have been merged.
The couring and main diagonal pointer is now manged by the patternm which means that they are calculated once only even if the preconditioner is deleted.



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