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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2470 - (show annotations)
Thu Jun 11 08:32:32 2009 UTC (10 years, 3 months ago) by gross
File MIME type: text/x-python
File size: 157691 byte(s)
SolverOption class added. It is not linked with anything.
1
2 ########################################################
3 #
4 # Copyright (c) 2003-2008 by University of Queensland
5 # Earth Systems Science Computational Center (ESSCC)
6 # http://www.uq.edu.au/esscc
7 #
8 # Primary Business: Queensland, Australia
9 # Licensed under the Open Software License version 3.0
10 # http://www.opensource.org/licenses/osl-3.0.php
11 #
12 ########################################################
13
14 __copyright__="""Copyright (c) 2003-2008 by University of Queensland
15 Earth Systems Science Computational Center (ESSCC)
16 http://www.uq.edu.au/esscc
17 Primary Business: Queensland, Australia"""
18 __license__="""Licensed under the Open Software License version 3.0
19 http://www.opensource.org/licenses/osl-3.0.php"""
20 __url__="https://launchpad.net/escript-finley"
21
22 """
23 The module provides an interface to define and solve linear partial
24 differential equations (PDEs) and Transport problems within L{escript}.
25 L{linearPDEs} does not provide any solver capabilities in itself but hands the
26 PDE over to the PDE solver library defined through the L{Domain<escript.Domain>}
27 of the PDE. The general interface is provided through the L{LinearPDE} class.
28 L{TransportProblem} provides an interface to initial value problems dominated
29 by its advective terms.
30
31 @var __author__: name of author
32 @var __copyright__: copyrights
33 @var __license__: licence agreement
34 @var __url__: url entry point on documentation
35 @var __version__: version
36 @var __date__: date of the version
37 """
38
39 import math
40 import escript
41 import util
42 import numpy
43
44 __author__="Lutz Gross, l.gross@uq.edu.au"
45
46
47 class SolverOptions(object):
48 """
49 this class defines the solver options for a linear or non-linear solver.
50
51 The option also supports the handling of diagnostic informations.
52
53 Typical usage is
54
55 opts=SolverOptions()
56 print opts
57 opts.resetDiagnostics()
58 u=solver(opts)
59 print "number of iteration steps: =",opts.getDiagnostics("num_iter")
60
61
62 @cvar DEFAULT: The default method used to solve the system of linear equations
63 @cvar DIRECT: The direct solver based on LDU factorization
64 @cvar CHOLEVSKY: The direct solver based on LDLt factorization (can only be applied for symmetric PDEs)
65 @cvar PCG: The preconditioned conjugate gradient method (can only be applied for symmetric PDEs)
66 @cvar CR: The conjugate residual method
67 @cvar CGS: The conjugate gradient square method
68 @cvar BICGSTAB: The stabilized Bi-Conjugate Gradient method
69 @cvar TFQMR: Transport Free Quasi Minimal Residual method
70 @cvar MINRES: Minimum residual method
71 @cvar SSOR: The symmetric over-relaxation method
72 @cvar ILU0: The incomplete LU factorization preconditioner with no fill-in
73 @cvar ILUT: The incomplete LU factorization preconditioner with fill-in
74 @cvar JACOBI: The Jacobi preconditioner
75 @cvar GMRES: The Gram-Schmidt minimum residual method
76 @cvar PRES20: Special GMRES with restart after 20 steps and truncation after 5 residuals
77 @cvar LUMPING: Matrix lumping
78 @cvar NO_REORDERING: No matrix reordering allowed
79 @cvar MINIMUM_FILL_IN: Reorder matrix to reduce fill-in during factorization
80 @cvar NESTED_DISSECTION: Reorder matrix to improve load balancing during factorization
81 @cvar PASO: PASO solver package
82 @cvar SCSL: SGI SCSL solver library
83 @cvar MKL: Intel's MKL solver library
84 @cvar UMFPACK: The UMFPACK library
85 @cvar TRILINOS: The TRILINOS parallel solver class library from Sandia National Labs
86 @cvar ITERATIVE: The default iterative solver
87 @cvar AMG: Algebraic Multi Grid
88 @cvar REC_ILU: recursive ILU0
89 @cvar RILU: relaxed ILU0
90 @cvar GAUSS_SEIDEL: Gauss-Seidel solver
91 @cvar DEFAULT_REORDERING: the reordering method recommended by the solver
92 @cvar SUPER_LU: the Super_LU solver package
93 @cvar PASTIX: the Pastix direct solver_package
94 @cvar YAIR_SHAPIRA_COARSENING: AMG coarsening method by Yair-Shapira
95 @cvar RUGE_STUEBEN_COARSENING: AMG coarsening method by Ruge and Stueben
96 @cvar AGGREGATION_COARSENING: AMG coarsening using (symmetric) aggregation
97 @cvar NO_PRECONDITIONER: no preconditioner is applied.
98 """
99 DEFAULT= 0
100 DIRECT= 1
101 CHOLEVSKY= 2
102 PCG= 3
103 CR= 4
104 CGS= 5
105 BICGSTAB= 6
106 SSOR= 7
107 ILU0= 8
108 ILUT= 9
109 JACOBI= 10
110 GMRES= 11
111 PRES20= 12
112 LUMPING= 13
113 NO_REORDERING= 17
114 MINIMUM_FILL_IN= 18
115 NESTED_DISSECTION= 19
116 MKL= 15
117 UMFPACK= 16
118 ITERATIVE= 20
119 PASO= 21
120 AMG= 22
121 REC_ILU = 23
122 TRILINOS = 24
123 NONLINEAR_GMRES = 25
124 TFQMR = 26
125 MINRES = 27
126 GAUSS_SEIDEL=28
127 RILU=29
128 DEFAULT_REORDERING=30
129 SUPER_LU=31
130 PASTIX=32
131 YAIR_SHAPIRA_COARSENING=33
132 RUGE_STUEBEN_COARSENING=34
133 AGGREGATION_COARSENING=35
134 NO_PRECONDITIONER=36
135 def __init__(self):
136 self.setLevelMax()
137 self.setCoarseningThreshold()
138 self.setNumSweeps()
139 self.setNumPreSweeps()
140 self.setNumPostSweeps()
141 self.setTolerance()
142 self.setAbsoluteTolerance()
143 self.setInnerTolerance()
144 self.setDropTolerance()
145 self.setDropStorage()
146 self.setIterMax()
147 self.setInnerIterMax()
148 self.setTruncation()
149 self.setRestart()
150 self.setSymmetry()
151 self.setVerbosity()
152 self.setInnerToleranceAdaption()
153 self.setAcceptanceConvergenceFailure()
154 self.setReordering()
155 self.setPackage()
156 self.setSolverMethod()
157 self.setPreconditioner()
158 self.setCoarsening()
159 self.setRelaxationFactor()
160 self.resetDiagnostics(all=True)
161
162 def __str__(self):
163 return self.getSummary()
164 def getSummary(self):
165 """
166 Returns a string reporting the current settings
167 """
168 out="Solver Package: %s"%(self.getName(self.getPackage()))
169 out+="\nVerbosity = %s"%self.isVerbose()
170 out+="\nAccept failed convergence = %s"%self.acceptConvergenceFailure()
171 out+="\nRelative tolerance = %e"%self.getTolerance()
172 out+="\nAbsolute tolerance = %e"%self.getAbsoluteTolerance()
173 out+="\nSymmetric problem = %s"%self.isSymmetric()
174 out+="\nMaximum number of iteration steps = %s"%self.getIterMax()
175 if self.getPackage() == self.PASO:
176 out+="\nSolver method = %s"%self.getName(self.getSolverMethod())
177 if self.getSolverMethod() == self.GMRES:
178 out+="\nTruncation = %s"%self.getTruncation()
179 out+="\nRestart = %s"%self.getRestart()
180 if self.getSolverMethod() == self.AMG:
181 out+="\nNumber of pre / post sweeps = %s / %s, %s"%(self.getNumPreSweeps(), self.getNumPostSweeps(), self.getNumSweeps())
182 out+="\nMaximum number of levels = %s"%self.LevelMax()
183 out+="\nCoarsening threshold = %e"%self.getCoarseningThreshold()
184 out+="\Coarsening method = %s"%self.getName(self.getCoarsening())
185 out+="\nPreconditioner = %s"%self.getName(self.getPreconditioner())
186 if self.getPreconditioner() == self.AMG:
187 out+="\nMaximum number of levels = %s"%self.LevelMax()
188 out+="\nCoarsening method = %s"%self.getName(self.getCoarsening())
189 out+="\nCoarsening threshold = %e"%self.getCoarseningThreshold()
190 out+="\nNumber of pre / post sweeps = %s / %s, %s"%(self.getNumPreSweeps(), self.getNumPostSweeps(), self.getNumSweeps())
191 if self.getPreconditioner() == self.GAUSS_SEIDEL:
192 out+="\nNumber of sweeps = %s"%self.getNumSweeps()
193 if self.getPreconditioner() == self.ILUT:
194 out+="\nDrop tolerance = %e"%self.getDropTolerance()
195 out+="\nStorage increase = %e"%self.getDropStorage()
196 if self.getPreconditioner() == self.RILU:
197 out+="\nRelaxation factor = %e"%self.getRelaxationFactor()
198 return out
199
200 def getName(self,key):
201 """
202 returns the name of a given key
203
204 @param key: a valid key
205 """
206 if key == self.DEFAULT: return "DEFAULT"
207 if key == self.DIRECT: return "DIRECT"
208 if key == self.CHOLEVSKY: return "CHOLEVSKY"
209 if key == self.PCG: return "PCG"
210 if key == self.CR: return "CR"
211 if key == self.CGS: return "CGS"
212 if key == self.BICGSTAB: return "BICGSTAB"
213 if key == self.SSOR: return "SSOR"
214 if key == self.ILU0: return "ILU0:"
215 if key == self.ILUT: return "ILUT"
216 if key == self.JACOBI: return "JACOBI"
217 if key == self.GMRES: return "GMRES"
218 if key == self.PRES20: return "PRES20"
219 if key == self.LUMPING: return "LUMPING"
220 if key == self.NO_REORDERING: return "NO_REORDERING"
221 if key == self.MINIMUM_FILL_IN: return "MINIMUM_FILL_IN"
222 if key == self.NESTED_DISSECTION: return "NESTED_DISSECTION"
223 if key == self.MKL: return "MKL"
224 if key == self.UMFPACK: return "UMFPACK"
225 if key == self.ITERATIVE: return "ITERATIVE"
226 if key == self.PASO: return "PASO"
227 if key == self.AMG: return "AMG"
228 if key == self.REC_ILU: return "REC_ILU"
229 if key == self.TRILINOS: return "TRILINOS"
230 if key == self.NONLINEAR_GMRES: return "NONLINEAR_GMRES"
231 if key == self.TFQMR: return "TFQMR"
232 if key == self.MINRES: return "MINRES"
233 if key == self.GAUSS_SEIDEL: return "GAUSS_SEIDEL"
234 if key == self.RILU: return "RILU"
235 if key == self.DEFAULT_REORDERING: return "DEFAULT_REORDERING"
236 if key == self.SUPER_LU: return "SUPER_LU"
237 if key == self.PASTIX: return "PASTIX"
238 if key == self.YAIR_SHAPIRA_COARSENING: return "YAIR_SHAPIRA_COARSENING"
239 if key == self.RUGE_STUEBEN_COARSENING: return "RUGE_STUEBEN_COARSENING"
240 if key == self.AGGREGATION_COARSENING: return "AGGREGATION_COARSENING"
241 if key == self.NO_PRECONDITIONER: return "NO_PRECONDITIONER"
242
243 def resetDiagnostics(self,all=False):
244 """
245 resets the diagnostics
246
247 @param all: if C{all} is C{True} all diagnostics including accumulative counters are reset.
248 @type all: C{bool}
249 """
250 self.__num_iter=None
251 self.__num_level=None
252 self.__num_inner_iter=None
253 self.__time=None
254 self.__set_up_time=None
255 self.__residual_norm=None
256 self.__converged=None
257 if all:
258 self.__cum_num_inner_iter=0
259 self.__cum_num_iter=0
260 self.__cum_time=0
261 self.__cum_set_up_time=0
262
263 def _updateDiagnostics(self, name, value):
264 """
265 Updates diagnostic information
266
267 @param name: name of diagnostic information
268 @type name: C{str} in the list "num_iter", "num_level", "num_inner_iter", "time", "set_up_time", "residual_norm", "converged".
269 @param vale: new value of the diagnostic information
270 @note: this function is used by a solver to report diagnostics informations.
271 """
272 if name == "num_iter":
273 self.__num_iter=int(value)
274 self.__cum_num_iter+=self.__num_iter
275 if name == "num_level":
276 self.__num_iter=int(value)
277 if name == "num_inner_iter":
278 self.__num_inner_iter=int(value)
279 self.__cum_num_inner_iter+=self.__num_inner_iter
280 if name == "time":
281 self.__time=float(value)
282 self.__cum_time+=self.__time
283 if name == "set_up_time":
284 self.__set_up_time=float(value)
285 self.__cum_set_up_time+=self.__set_up_time
286 if name == "residual_norm":
287 self.__residual_norm=float(value)
288 if name == "converged":
289 self.__converged = (value == True)
290 def getDiagnostics(self, name):
291 """
292 Returns the diagnostic information C{name}
293
294 @param name: name of diagnostic information where
295 - "num_iter": the number of iteration steps
296 - "cum_num_iter": the cumulative number of iteration steps
297 - "num_level": the number of level in multi level solver
298 - "num_inner_iter": the number of inner iteration steps
299 - "cum_num_inner_iter": the cumulative number of inner iteration steps
300 - "time": execution time
301 - "cum_time": cumulative execution time
302 - "set_up_time": time to set up of the solver, typically this includes factorization and reordering
303 - "cum_set_up_time": cumulative time to set up of the solver
304 - "residual_norm": norm of the final residual
305 - "converged": return self.__converged
306 @type name: C{str} in the list "num_iter", "num_level", "num_inner_iter", "time", "set_up_time", "residual_norm", "converged".
307 @return: requested value. C{None} is returned if the value is undefined.
308 @note: If the solver has thrown an exception diagnostic values have an undefined status.
309 """
310 if name == "num_iter": return self.__num_iter
311 elif name == "cum_num_iter": return self.__cum_num_iter
312 elif name == "num_level": return self.__num_level
313 elif name == "num_inner_iter": return self.__num_inner_iter
314 elif name == "cum_num_inner_iter": return self.__cum_num_inner_iter
315 elif name == "time": return self.__time
316 elif name == "cum_time": return self.__cum_time
317 elif name == "set_up_time": return self.__set_up_time
318 elif name == "cum_set_up_time": return self.__cum_set_up_time
319 elif name == "residual_norm": return self.__residual_norm
320 elif name == "converged": return self.__converged
321 else:
322 raise ValueError,"unknown diagnostic item %s"%name
323 def hasConverged(self):
324 """
325 Returns C{True} if the last solver call has been finalized successfully.
326 @note: if an exception has been thrown by the solver the status of this flag is undefined.
327 """
328 return self.getDiagnostics("converged")
329 def setCoarsening(self,method=0):
330 """
331 Sets the key of the coarsening method to be applied in AMG.
332
333 @param method: selects the coarsening method .
334 @type method: in L{SolverOptions.DEFAULT}, L{SolverOptions.YAIR_SHAPIRA_COARSENING},
335 L{SolverOptions.RUGE_STUEBEN_COARSENING}, L{SolverOptions.AGGREGATION_COARSENING}
336 """
337 if not method in [self.DEFAULT, self.YAIR_SHAPIRA_COARSENING, self.RUGE_STUEBEN_COARSENING, self.AGGREGATION_COARSENING]:
338 raise ValueError,"unknown coarsening method %s"%method
339 self.__coarsening=method
340 def getCoarsening(self):
341 """
342 Returns the key of the coarsening algorithm to be applied AMG.
343
344 @rtype: in the list L{SolverOptions.DEFAULT}, L{SolverOptions.YAIR_SHAPIRA_COARSENING},
345 L{SolverOptions.RUGE_STUEBEN_COARSENING}, L{SolverOptions.AGGREGATION_COARSENING}
346 """
347 return self.__coarsening
348 def setPreconditioner(self, preconditioner=10):
349 """
350 Sets the preconditioner to be used.
351
352 @param preconditioner: key of the preconditioner to be used.
353 @type preconditioner: in L{SolverOptions.SSOR}, L{SolverOptions.ILU0}, L{SolverOptions.ILUT}, L{SolverOptions.JACOBI},
354 L{SolverOptions.AMG}, L{SolverOptions.REC_ILU}, L{SolverOptions.GAUSS_SEIDEL}, L{SolverOptions.RILU},
355 L{SolverOptions.NO_PRECONDITIONER}
356 @note: Not all packages support all preconditioner. It can be assumed that a package makes a reasonable choice if it encounters
357 an unknown preconditioner.
358 """
359 if not preconditioner in [ SolverOptions.SSOR, SolverOptions.ILU0, SolverOptions.ILUT, SolverOptions.JACOBI,
360 SolverOptions.AMG, SolverOptions.REC_ILU, SolverOptions.GAUSS_SEIDEL, SolverOptions.RILU,
361 SolverOptions.NO_PRECONDITIONER] :
362 raise ValueError,"unknown preconditioner %s"%preconditioner
363 self.__preconditioner=preconditioner
364 def getPreconditioner(self):
365 """
366 Returns key of the preconditioner to be used.
367
368 @rtype: in the list L{SolverOptions.SSOR}, L{SolverOptions.ILU0}, L{SolverOptions.ILUT}, L{SolverOptions.JACOBI},
369 L{SolverOptions.AMG}, L{SolverOptions.REC_ILU}, L{SolverOptions.GAUSS_SEIDEL}, L{SolverOptions.RILU},
370 L{SolverOptions.NO_PRECONDITIONER}
371 """
372 return self.__preconditioner
373 def setSolverMethod(self, method=0):
374 """
375 Sets the solver method to be used. Use C{method}=C{SolverOptions.DIRECT} to indicate that a direct rather than an iterative
376 solver should be used and Use C{method}=C{SolverOptions.ITERATIVE} to indicate that an iterative rather than a direct
377 solver should be used.
378
379 @param method: key of the solver method to be used.
380 @type method: in L{SolverOptions.DEFAULT}, L{SolverOptions.DIRECT}, L{SolverOptions.CHOLEVSKY}, L{SolverOptions.PCG},
381 L{SolverOptions.CR}, L{SolverOptions.CGS}, L{SolverOptions.BICGSTAB}, L{SolverOptions.SSOR},
382 L{SolverOptions.GMRES}, L{SolverOptions.PRES20}, L{SolverOptions.LUMPING}, L{SolverOptions.ITERATIVE},
383 L{SolverOptions.AMG}, L{SolverOptions.NONLINEAR_GMRES}, L{SolverOptions.TFQMR}, L{SolverOptions.MINRES},
384 L{SolverOptions.GAUSS_SEIDEL}
385 @note: Not all packages support all solvers. It can be assumed that a package makes a reasonable choice if it encounters
386 an unknown solver method.
387 """
388 if not method in [ SolverOptions.DEFAULT, SolverOptions.DIRECT, SolverOptions.CHOLEVSKY, SolverOptions.PCG,
389 SolverOptions.CR, SolverOptions.CGS, SolverOptions.BICGSTAB, SolverOptions.SSOR,
390 SolverOptions.GMRES, SolverOptions.PRES20, SolverOptions.LUMPING, SolverOptions.ITERATIVE, SolverOptions.AMG,
391 SolverOptions.NONLINEAR_GMRES, SolverOptions.TFQMR, SolverOptions.MINRES, SolverOptions.GAUSS_SEIDEL]:
392 raise ValueError,"unknown solver method %s"%method
393 self.__method=method
394 def getSolverMethod(self):
395 """
396 Returns key of the solver method to be used.
397
398 @rtype: in the list L{SolverOptions.DEFAULT}, L{SolverOptions.DIRECT}, L{SolverOptions.CHOLEVSKY}, L{SolverOptions.PCG},
399 L{SolverOptions.CR}, L{SolverOptions.CGS}, L{SolverOptions.BICGSTAB}, L{SolverOptions.SSOR},
400 L{SolverOptions.GMRES}, L{SolverOptions.PRES20}, L{SolverOptions.LUMPING}, L{SolverOptions.ITERATIVE},
401 L{SolverOptions.AMG}, L{SolverOptions.NONLINEAR_GMRES}, L{SolverOptions.TFQMR}, L{SolverOptions.MINRES},
402 L{SolverOptions.GAUSS_SEIDEL}
403 """
404 return self.__method
405
406 def setPackage(self, package=0):
407 """
408 Sets the solver package to be used as a solver.
409
410 @param package: key of the solver package to be used.
411 @type package: in L{SolverOptions.DEFAULT}, L{SolverOptions.PASO}, L{SolverOptions.SUPER_LU}, L{SolverOptions.PASTIX}, L{SolverOptions.MKL}, L{SolverOptions.UMFPACK}, L{SolverOptions.TRILINOS}
412 @note: Not all packages are support on all implementation. An exception may be thrown on some platforms if a particular is requested.
413 """
414 if not package in [SolverOptions.DEFAULT, SolverOptions.PASO, SolverOptions.SUPER_LU, SolverOptions.PASTIX, SolverOptions.MKL, SolverOptions.UMFPACK, SolverOptions.TRILINOS]:
415 raise ValueError,"unknown solver package %s"%package
416 self.__package=package
417 def getPackage(self):
418 """
419 Returns the solver package key
420
421 @rtype: in the list L{SolverOptions.DEFAULT}, L{SolverOptions.PASO}, L{SolverOptions.SUPER_LU}, L{SolverOptions.PASTIX}, L{SolverOptions.MKL}, L{SolverOptions.UMFPACK}, L{SolverOptions.TRILINOS}
422 """
423 return self.__package
424 def setReordering(self,ordering=30):
425 """
426 Sets the key of the reordering method to be applied if supported by the solver. Some direct solvers support reordering
427 to optimize compute time and storage use during elimination.
428
429 @param ordering: selects the reordering strategy.
430 @type ordering: in L{SolverOptions.NO_REORDERING}, L{SolverOptions.NO_REORDERING},
431 L{SolverOptions.NO_REORDERING}, L{SolverOptions.DEFAULT_REORDERING}
432 """
433 if not ordering in [self.NO_REORDERING, self.MINIMUM_FILL_IN, self.NESTED_DISSECTION, self.DEFAULT_REORDERING]:
434 raise ValueError,"unknown reordering strategy %s"%ordering
435 self.__reordering=ordering
436 def getReordering(self):
437 """
438 Returns the key of the reordering method to be applied if supported by the solver.
439
440 @rtype: in the list L{SolverOptions.NO_REORDERING}, L{SolverOptions.NO_REORDERING},
441 L{SolverOptions.NO_REORDERING}, L{SolverOptions.DEFAULT_REORDERING}
442 """
443 return self.__reordering
444 def setRestart(self,restart=None):
445 """
446 Sets the number of iterations steps after which GMRES is performing a restart.
447
448 @param restart: number of iteration steps after which to perform a restart. If equal to C{None} no
449 restart is performed.
450 @type restart: C{int} or C{None}
451 """
452 if restart == None:
453 self.__restart=restart
454 else:
455 restart=int(restart)
456 if restart<1:
457 raise ValueError,"restart must be positive."
458 self.__restart=restart
459 def getRestart(self):
460 """
461 Returns the number of iterations steps after which GMRES is performing a restart.
462 If C{None} is returned no restart is performed.
463
464 @rtype: C{int} or C{None}
465 """
466 if self.__restart < 0:
467 return None
468 else:
469 return self.__restart
470 def setTruncation(self,truncation=20):
471 """
472 Sets the number of residuals in GMRES to be stored for orthogonalization. The more residuals are stored
473 the faster GMRES converged but
474
475 @param truncation: truncation
476 @type truncation: C{int}
477 """
478 truncation=int(truncation)
479 if truncation<1:
480 raise ValueError,"truncation must be positive."
481 self.__truncation=truncation
482 def getTruncation(self):
483 """
484 Returns the number of residuals in GMRES to be stored for orthogonalization
485
486 @rtype: C{int}
487 """
488 return self.__truncation
489 def setInnerIterMax(self,iter_max=10):
490 """
491 Sets the maximum number of iteration steps for the inner iteration.
492
493 @param iter_max: maximum number of inner iterations
494 @type iter_max: C{int}
495 """
496 iter_max=int(iter_max)
497 if iter_max<1:
498 raise ValueError,"maximum number of inner iteration must be positive."
499 self.__inner_iter_max=iter_max
500 def getInnerIterMax(self):
501 """
502 Returns maximum number of inner iteration steps
503
504 @rtype: C{int}
505 """
506 return self.__inner_iter_max
507 def setIterMax(self,iter_max=10000):
508 """
509 Sets the maximum number of iteration steps
510
511 @param iter_max: maximum number of iteration steps
512 @type iter_max: C{int}
513 """
514 iter_max=int(iter_max)
515 if iter_max<1:
516 raise ValueError,"maximum number of iteration steps must be positive."
517 self.__iter_max=iter_max
518 def getIterMax(self):
519 """
520 Returns maximum number of iteration steps
521
522 @rtype: C{int}
523 """
524 return self.__iter_max
525 def setLevelMax(self,level_max=10):
526 """
527 Sets the maximum number of coarsening levels to be used in an algebraic multi level solver or preconditioner
528
529 @param level_max: maximum number of levels
530 @type level_max: C{int}
531 """
532 level_max=int(level_max)
533 if level_max<0:
534 raise ValueError,"maximum number of coarsening levels must be non-negative."
535 self.__level_max=level_max
536 def getLevelMax(self):
537 """
538 Returns the maximum number of coarsening levels to be used in an algebraic multi level solver or preconditioner
539
540 @rtype: C{int}
541 """
542 return self.__level_max
543 def setCoarseningThreshold(self,theta=0.05):
544 """
545 Sets the threshold for coarsening in the algebraic multi level solver or preconditioner
546
547 @param theta: threshold for coarsening
548 @type theta: positive C{float}
549 """
550 theta=float(theta)
551 if theta<0 or theta>1:
552 raise ValueError,"threshold must be non-negative and less or equal 1."
553 self.__coarsening_threshold=theta
554 def getCoarseningThreshold(self):
555 """
556 Returns the threshold for coarsening in the algebraic multi level solver or preconditioner
557
558 @rtype: C{float}
559 """
560 return self.__coarsening_threshold
561 def setNumSweeps(self,sweeps=2):
562 """
563 Sets the number of sweeps in a Jacobi or Gauss-Seidel/SOR preconditioner.
564
565 @param sweeps: number of sweeps
566 @type theta: positive C{int}
567 """
568 sweeps=int(sweeps)
569 if sweeps<1:
570 raise ValueError,"number of sweeps must be positive."
571 self.__sweeps=sweeps
572 def getNumSweeps(self):
573 """
574 Returns the number of sweeps in a Jacobi or Gauss-Seidel/SOR preconditioner.
575
576 @rtype: C{int}
577 """
578 return self.__sweeps
579 def setNumPreSweeps(self,sweeps=2):
580 """
581 Sets the number of sweeps in the pre-smoothing step of a multi level solver or preconditioner
582
583 @param sweeps: number of sweeps
584 @type theta: positive C{int}
585 """
586 sweeps=int(sweeps)
587 if sweeps<1:
588 raise ValueError,"number of sweeps must be positive."
589 self.__pre_sweeps=sweeps
590 def getNumPreSweeps(self):
591 """
592 Returns he number of sweeps in the pre-smoothing step of a multi level solver or preconditioner
593
594 @rtype: C{int}
595 """
596 return self.__pre_sweeps
597 def setNumPostSweeps(self,sweeps=2):
598 """
599 Sets the number of sweeps in the post-smoothing step of a multi level solver or preconditioner
600
601 @param sweeps: number of sweeps
602 @type theta: positive C{int}
603 """
604 sweeps=int(sweeps)
605 if sweeps<1:
606 raise ValueError,"number of sweeps must be positive."
607 self.__post_sweeps=sweeps
608 def getNumPostSweeps(self):
609 """
610 Returns he number of sweeps in the post-smoothing step of a multi level solver or preconditioner
611
612 @rtype: C{int}
613 """
614 return self.__post_sweeps
615
616 def setTolerance(self,rtol=1.e-8):
617 """
618 Sets the relative tolerance for the solver
619
620 @param rtol: relative tolerance
621 @type rtol: non-negative C{float}
622 """
623 rtol=float(rtol)
624 if rtol<0 or rtol>1:
625 raise ValueError,"tolerance must be non-negative and less or equal 1."
626 self.__tolerance=rtol
627 def getTolerance(self):
628 """
629 Returns the relative tolerance for the solver
630
631 @rtype: C{float}
632 """
633 return self.__tolerance
634 def setAbsoluteTolerance(self,atol=0.):
635 """
636 Sets the absolute tolerance for the solver
637
638 @param atol: absolute tolerance
639 @type atol: non-negative C{float}
640 """
641 atol=float(atol)
642 if atol<0:
643 raise ValueError,"tolerance must be non-negative."
644 self.__absolute_tolerance=atol
645 def getAbsoluteTolerance(self):
646 """
647 Returns the absolute tolerance for the solver
648
649 @rtype: C{float}
650 """
651 return self.__absolute_tolerance
652
653 def setInnerTolerance(self,rtol=0.9):
654 """
655 Sets the relative tolerance for an inner iteration scheme for instance
656 on the coarsest level in a multi-level scheme.
657
658 @param rtol: inner relative tolerance
659 @type rtol: positive C{float}
660 """
661 rtol=float(rtol)
662 if rtol<=0 or rtol>1:
663 raise ValueError,"tolerance must be positive and less or equal 1."
664 self.__inner_tolerance=rtol
665 def getInnerTolerance(self):
666 """
667 Returns the relative tolerance for an inner iteration scheme
668
669 @rtype: C{float}
670 """
671 return self.__inner_tolerance
672 def setDropTolerance(self,drop_tol=0.01):
673 """
674 Sets the relative drop tolerance in ILUT
675
676 @param drop_tol: drop tolerance
677 @type drop_tol: positive C{float}
678 """
679 drop_tol=float(drop_tol)
680 if drop_tol<=0 or drop_tol>1:
681 raise ValueError,"drop tolerance must be positive and less or equal 1."
682 self.__drop_tolerance=drop_tol
683 def getDropTolerance(self):
684 """
685 Returns the relative drop tolerance in ILUT
686
687 @rtype: C{float}
688 """
689 return self.__drop_tolerance
690 def setDropStorage(self,storage=2.):
691 """
692 Sets the maximum allowed increase in storage for ILUT. C{storage}=2 would mean that
693 a doubling of the storage needed for the coefficient matrix is allowed in the ILUT factorization.
694
695 @param storage: allowed storage increase
696 @type storage: C{float}
697 """
698 storage=float(storage)
699 if storage<1:
700 raise ValueError,"allowed storage increase must be greater or equal to 1."
701 self.__drop_storage=storage
702 def getDropStorage(self):
703
704 """
705 Returns the maximum allowed increase in storage for ILUT
706
707 @rtype: C{float}
708 """
709 return self.__drop_storage
710 def setRelaxationFactor(self,factor=0.3):
711 """
712 Sets the relaxation factor used to add dropped elements in RILU to the main diagonal.
713
714 @param factor: relaxation factor
715 @type factor: C{float}
716 @note: RILU with a relaxation factor 0 is identical to ILU0
717 """
718 factor=float(factor)
719 if factor<0:
720 raise ValueError,"relaxation factor must be non-negative."
721 self.__relaxation=factor
722 def getRelaxationFactor(self):
723
724 """
725 Returns the relaxation factor used to add dropped elements in RILU to the main diagonal.
726
727 @rtype: C{float}
728 """
729 return self.__relaxation
730 def isSymmetric(self):
731 """
732 Checks if symmetry of the coefficient matrix is indicated.
733
734 @return: True if a symmetric PDE is indicated, False otherwise
735 @rtype: C{bool}
736 """
737 return self.__symmetric
738 def setSymmetryOn(self):
739 """
740 Sets the symmetry flag to indicate that the coefficient matrix is symmetric.
741 """
742 self.__symmetric=True
743 def setSymmetryOff(self):
744 """
745 Clears the symmetry flag for the coefficient matrix.
746 """
747 self.__symmetric=False
748 def setSymmetry(self,flag=False):
749 """
750 Sets the symmetry flag for the coefficient matrix to C{flag}.
751
752 @param flag: If True, the symmetry flag is set otherwise reset.
753 @type flag: C{bool}
754 """
755 if flag:
756 self.setSymmetryOn()
757 else:
758 self.setSymmetryOff()
759 def isVerbose(self):
760 """
761 Returns C{True} if the solver is expected to be verbose.
762
763 @return: True if verbosity of switched on.
764 @rtype: C{bool}
765 """
766 return self.__verbose
767
768 def setVerbosityOn(self):
769 """
770 Switches the verbosity of the solver on.
771 """
772 self.__verbose=True
773 def setVerbosityOff(self):
774 """
775 Switches the verbosity of the solver off.
776 """
777 self.__verbose=False
778 def setVerbosity(self,verbose=True):
779 """
780 Sets the verbosity flag for the solver to C{flag}.
781
782 @param flag: If C{True}, the verbosity of the solver is switched on.
783 @type flag: C{bool}
784 """
785 if verbose:
786 self.setVerbosityOn()
787 else:
788 self.setVerbosityOff()
789
790 self.__adapt_inner_tolerance=True
791 def adaptInnerTolerance(self):
792 """
793 Returns C{True} if the tolerance of the inner solver is selected automatically.
794 Otherwise the inner tolerance set by L{setInnerTolerance} is used.
795
796 @return: C{True} if inner tolerance adaption is chosen.
797 @rtype: C{bool}
798 """
799 return self.__adapt_inner_tolerance
800
801 def setInnerToleranceAdaptionOn(self):
802 """
803 Switches the automatic selection of inner tolerance on
804 """
805 self.__adapt_inner_tolerance=True
806 def setInnerToleranceAdaptionOff(self):
807 """
808 Switches the automatic selection of inner tolerance off.
809 """
810 self.__adapt_inner_tolerance=False
811 def setInnerToleranceAdaption(self,adapt=True):
812 """
813 Sets a flag to indicate automatic selection of the inner tolerance.
814
815 @param adapt: If C{True}, the inner tolerance is selected automatically.
816 @type adapt: C{bool}
817 """
818 if adapt:
819 self.setInnerToleranceAdaptionOn()
820 else:
821 self.setInnerToleranceAdaptionOff()
822
823 def acceptConvergenceFailure(self):
824 """
825 Returns C{True} if a failure to meet the stopping criteria within the
826 given number of iteration steps is not raising in exception. This is useful
827 if a solver is used in a non-linear context where the non-linear solver can
828 continue even if the returned the solution does not necessarily meet the
829 stopping criteria. One can use the L{hasConverged} method to check if the
830 last call to the solver was successful.
831
832 @return: C{True} if a failure to achieve convergence is accepted.
833 @rtype: C{bool}
834 """
835 return self.__accept_convergence_failure
836
837 def setAcceptanceConvergenceFailureOn(self):
838 """
839 Switches the acceptance of a failure of convergence on
840 """
841 self.__accept_convergence_failure=True
842 def setAcceptanceConvergenceFailureOff(self):
843 """
844 Switches the acceptance of a failure of convergence off.
845 """
846 self.__accept_convergence_failure=False
847 def setAcceptanceConvergenceFailure(self,accept=False):
848 """
849 Sets a flag to indicate the acceptance of a failure of convergence.
850
851 @param accept: If C{True}, any failure to achieve convergence is accepted.
852 @type accept: C{bool}
853 """
854 if accept:
855 self.setAcceptanceConvergenceFailureOn()
856 else:
857 self.setAcceptanceConvergenceFailureOff()
858
859 class IllegalCoefficient(ValueError):
860 """
861 Exception that is raised if an illegal coefficient of the general or
862 particular PDE is requested.
863 """
864 pass
865
866 class IllegalCoefficientValue(ValueError):
867 """
868 Exception that is raised if an incorrect value for a coefficient is used.
869 """
870 pass
871
872 class IllegalCoefficientFunctionSpace(ValueError):
873 """
874 Exception that is raised if an incorrect function space for a coefficient
875 is used.
876 """
877
878 class UndefinedPDEError(ValueError):
879 """
880 Exception that is raised if a PDE is not fully defined yet.
881 """
882 pass
883
884 class PDECoef(object):
885 """
886 A class for describing a PDE coefficient.
887
888 @cvar INTERIOR: indicator that coefficient is defined on the interior of
889 the PDE domain
890 @cvar BOUNDARY: indicator that coefficient is defined on the boundary of
891 the PDE domain
892 @cvar CONTACT: indicator that coefficient is defined on the contact region
893 within the PDE domain
894 @cvar INTERIOR_REDUCED: indicator that coefficient is defined on the
895 interior of the PDE domain using a reduced
896 integration order
897 @cvar BOUNDARY_REDUCED: indicator that coefficient is defined on the
898 boundary of the PDE domain using a reduced
899 integration order
900 @cvar CONTACT_REDUCED: indicator that coefficient is defined on the contact
901 region within the PDE domain using a reduced
902 integration order
903 @cvar SOLUTION: indicator that coefficient is defined through a solution of
904 the PDE
905 @cvar REDUCED: indicator that coefficient is defined through a reduced
906 solution of the PDE
907 @cvar BY_EQUATION: indicator that the dimension of the coefficient shape is
908 defined by the number of PDE equations
909 @cvar BY_SOLUTION: indicator that the dimension of the coefficient shape is
910 defined by the number of PDE solutions
911 @cvar BY_DIM: indicator that the dimension of the coefficient shape is
912 defined by the spatial dimension
913 @cvar OPERATOR: indicator that the the coefficient alters the operator of
914 the PDE
915 @cvar RIGHTHANDSIDE: indicator that the the coefficient alters the right
916 hand side of the PDE
917 @cvar BOTH: indicator that the the coefficient alters the operator as well
918 as the right hand side of the PDE
919
920 """
921 INTERIOR=0
922 BOUNDARY=1
923 CONTACT=2
924 SOLUTION=3
925 REDUCED=4
926 BY_EQUATION=5
927 BY_SOLUTION=6
928 BY_DIM=7
929 OPERATOR=10
930 RIGHTHANDSIDE=11
931 BOTH=12
932 INTERIOR_REDUCED=13
933 BOUNDARY_REDUCED=14
934 CONTACT_REDUCED=15
935
936 def __init__(self, where, pattern, altering):
937 """
938 Initialises a PDE coefficient type.
939
940 @param where: describes where the coefficient lives
941 @type where: one of L{INTERIOR}, L{BOUNDARY}, L{CONTACT}, L{SOLUTION},
942 L{REDUCED}, L{INTERIOR_REDUCED}, L{BOUNDARY_REDUCED},
943 L{CONTACT_REDUCED}
944 @param pattern: describes the shape of the coefficient and how the shape
945 is built for a given spatial dimension and numbers of
946 equations and solutions in then PDE. For instance,
947 (L{BY_EQUATION},L{BY_SOLUTION},L{BY_DIM}) describes a
948 rank 3 coefficient which is instantiated as shape (3,2,2)
949 in case of three equations and two solution components
950 on a 2-dimensional domain. In the case of single equation
951 and a single solution component the shape components
952 marked by L{BY_EQUATION} or L{BY_SOLUTION} are dropped.
953 In this case the example would be read as (2,).
954 @type pattern: C{tuple} of L{BY_EQUATION}, L{BY_SOLUTION}, L{BY_DIM}
955 @param altering: indicates what part of the PDE is altered if the
956 coefficient is altered
957 @type altering: one of L{OPERATOR}, L{RIGHTHANDSIDE}, L{BOTH}
958 """
959 super(PDECoef, self).__init__()
960 self.what=where
961 self.pattern=pattern
962 self.altering=altering
963 self.resetValue()
964
965 def resetValue(self):
966 """
967 Resets the coefficient value to the default.
968 """
969 self.value=escript.Data()
970
971 def getFunctionSpace(self,domain,reducedEquationOrder=False,reducedSolutionOrder=False):
972 """
973 Returns the L{FunctionSpace<escript.FunctionSpace>} of the coefficient.
974
975 @param domain: domain on which the PDE uses the coefficient
976 @type domain: L{Domain<escript.Domain>}
977 @param reducedEquationOrder: True to indicate that reduced order is used
978 to represent the equation
979 @type reducedEquationOrder: C{bool}
980 @param reducedSolutionOrder: True to indicate that reduced order is used
981 to represent the solution
982 @type reducedSolutionOrder: C{bool}
983 @return: L{FunctionSpace<escript.FunctionSpace>} of the coefficient
984 @rtype: L{FunctionSpace<escript.FunctionSpace>}
985 """
986 if self.what==self.INTERIOR:
987 return escript.Function(domain)
988 elif self.what==self.INTERIOR_REDUCED:
989 return escript.ReducedFunction(domain)
990 elif self.what==self.BOUNDARY:
991 return escript.FunctionOnBoundary(domain)
992 elif self.what==self.BOUNDARY_REDUCED:
993 return escript.ReducedFunctionOnBoundary(domain)
994 elif self.what==self.CONTACT:
995 return escript.FunctionOnContactZero(domain)
996 elif self.what==self.CONTACT_REDUCED:
997 return escript.ReducedFunctionOnContactZero(domain)
998 elif self.what==self.SOLUTION:
999 if reducedEquationOrder and reducedSolutionOrder:
1000 return escript.ReducedSolution(domain)
1001 else:
1002 return escript.Solution(domain)
1003 elif self.what==self.REDUCED:
1004 return escript.ReducedSolution(domain)
1005
1006 def getValue(self):
1007 """
1008 Returns the value of the coefficient.
1009
1010 @return: value of the coefficient
1011 @rtype: L{Data<escript.Data>}
1012 """
1013 return self.value
1014
1015 def setValue(self,domain,numEquations=1,numSolutions=1,reducedEquationOrder=False,reducedSolutionOrder=False,newValue=None):
1016 """
1017 Sets the value of the coefficient to a new value.
1018
1019 @param domain: domain on which the PDE uses the coefficient
1020 @type domain: L{Domain<escript.Domain>}
1021 @param numEquations: number of equations of the PDE
1022 @type numEquations: C{int}
1023 @param numSolutions: number of components of the PDE solution
1024 @type numSolutions: C{int}
1025 @param reducedEquationOrder: True to indicate that reduced order is used
1026 to represent the equation
1027 @type reducedEquationOrder: C{bool}
1028 @param reducedSolutionOrder: True to indicate that reduced order is used
1029 to represent the solution
1030 @type reducedSolutionOrder: C{bool}
1031 @param newValue: number of components of the PDE solution
1032 @type newValue: any object that can be converted into a
1033 L{Data<escript.Data>} object with the appropriate shape
1034 and L{FunctionSpace<escript.FunctionSpace>}
1035 @raise IllegalCoefficientValue: if the shape of the assigned value does
1036 not match the shape of the coefficient
1037 @raise IllegalCoefficientFunctionSpace: if unable to interpolate value
1038 to appropriate function space
1039 """
1040 if newValue==None:
1041 newValue=escript.Data()
1042 elif isinstance(newValue,escript.Data):
1043 if not newValue.isEmpty():
1044 if not newValue.getFunctionSpace() == self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder):
1045 try:
1046 newValue=escript.Data(newValue,self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder))
1047 except:
1048 raise IllegalCoefficientFunctionSpace,"Unable to interpolate coefficient to function space %s"%self.getFunctionSpace(domain)
1049 else:
1050 newValue=escript.Data(newValue,self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder))
1051 if not newValue.isEmpty():
1052 if not self.getShape(domain,numEquations,numSolutions)==newValue.getShape():
1053 raise IllegalCoefficientValue,"Expected shape of coefficient is %s but actual shape is %s."%(self.getShape(domain,numEquations,numSolutions),newValue.getShape())
1054 self.value=newValue
1055
1056 def isAlteringOperator(self):
1057 """
1058 Checks if the coefficient alters the operator of the PDE.
1059
1060 @return: True if the operator of the PDE is changed when the
1061 coefficient is changed
1062 @rtype: C{bool}
1063 """
1064 if self.altering==self.OPERATOR or self.altering==self.BOTH:
1065 return not None
1066 else:
1067 return None
1068
1069 def isAlteringRightHandSide(self):
1070 """
1071 Checks if the coefficient alters the right hand side of the PDE.
1072
1073 @rtype: C{bool}
1074 @return: True if the right hand side of the PDE is changed when the
1075 coefficient is changed, C{None} otherwise.
1076 """
1077 if self.altering==self.RIGHTHANDSIDE or self.altering==self.BOTH:
1078 return not None
1079 else:
1080 return None
1081
1082 def estimateNumEquationsAndNumSolutions(self,domain,shape=()):
1083 """
1084 Tries to estimate the number of equations and number of solutions if
1085 the coefficient has the given shape.
1086
1087 @param domain: domain on which the PDE uses the coefficient
1088 @type domain: L{Domain<escript.Domain>}
1089 @param shape: suggested shape of the coefficient
1090 @type shape: C{tuple} of C{int} values
1091 @return: the number of equations and number of solutions of the PDE if
1092 the coefficient has given shape. If no appropriate numbers
1093 could be identified, C{None} is returned
1094 @rtype: C{tuple} of two C{int} values or C{None}
1095 """
1096 dim=domain.getDim()
1097 if len(shape)>0:
1098 num=max(shape)+1
1099 else:
1100 num=1
1101 search=[]
1102 if self.definesNumEquation() and self.definesNumSolutions():
1103 for u in range(num):
1104 for e in range(num):
1105 search.append((e,u))
1106 search.sort(self.__CompTuple2)
1107 for item in search:
1108 s=self.getShape(domain,item[0],item[1])
1109 if len(s)==0 and len(shape)==0:
1110 return (1,1)
1111 else:
1112 if s==shape: return item
1113 elif self.definesNumEquation():
1114 for e in range(num,0,-1):
1115 s=self.getShape(domain,e,0)
1116 if len(s)==0 and len(shape)==0:
1117 return (1,None)
1118 else:
1119 if s==shape: return (e,None)
1120
1121 elif self.definesNumSolutions():
1122 for u in range(num,0,-1):
1123 s=self.getShape(domain,0,u)
1124 if len(s)==0 and len(shape)==0:
1125 return (None,1)
1126 else:
1127 if s==shape: return (None,u)
1128 return None
1129
1130 def definesNumSolutions(self):
1131 """
1132 Checks if the coefficient allows to estimate the number of solution
1133 components.
1134
1135 @return: True if the coefficient allows an estimate of the number of
1136 solution components, False otherwise
1137 @rtype: C{bool}
1138 """
1139 for i in self.pattern:
1140 if i==self.BY_SOLUTION: return True
1141 return False
1142
1143 def definesNumEquation(self):
1144 """
1145 Checks if the coefficient allows to estimate the number of equations.
1146
1147 @return: True if the coefficient allows an estimate of the number of
1148 equations, False otherwise
1149 @rtype: C{bool}
1150 """
1151 for i in self.pattern:
1152 if i==self.BY_EQUATION: return True
1153 return False
1154
1155 def __CompTuple2(self,t1,t2):
1156 """
1157 Compares two tuples of possible number of equations and number of
1158 solutions.
1159
1160 @param t1: the first tuple
1161 @param t2: the second tuple
1162 @return: 0, 1, or -1
1163 """
1164
1165 dif=t1[0]+t1[1]-(t2[0]+t2[1])
1166 if dif<0: return 1
1167 elif dif>0: return -1
1168 else: return 0
1169
1170 def getShape(self,domain,numEquations=1,numSolutions=1):
1171 """
1172 Builds the required shape of the coefficient.
1173
1174 @param domain: domain on which the PDE uses the coefficient
1175 @type domain: L{Domain<escript.Domain>}
1176 @param numEquations: number of equations of the PDE
1177 @type numEquations: C{int}
1178 @param numSolutions: number of components of the PDE solution
1179 @type numSolutions: C{int}
1180 @return: shape of the coefficient
1181 @rtype: C{tuple} of C{int} values
1182 """
1183 dim=domain.getDim()
1184 s=()
1185 for i in self.pattern:
1186 if i==self.BY_EQUATION:
1187 if numEquations>1: s=s+(numEquations,)
1188 elif i==self.BY_SOLUTION:
1189 if numSolutions>1: s=s+(numSolutions,)
1190 else:
1191 s=s+(dim,)
1192 return s
1193
1194 #====================================================================================================================
1195
1196 class LinearProblem(object):
1197 """
1198 This is the base class to define a general linear PDE-type problem for
1199 for an unknown function M{u} on a given domain defined through a
1200 L{Domain<escript.Domain>} object. The problem can be given as a single
1201 equation or as a system of equations.
1202
1203 The class assumes that some sort of assembling process is required to form
1204 a problem of the form
1205
1206 M{L u=f}
1207
1208 where M{L} is an operator and M{f} is the right hand side. This operator
1209 problem will be solved to get the unknown M{u}.
1210
1211 @cvar DEFAULT: The default method used to solve the system of linear
1212 equations
1213 @cvar DIRECT: The direct solver based on LDU factorization
1214 @cvar CHOLEVSKY: The direct solver based on LDLt factorization (can only be
1215 applied for symmetric PDEs)
1216 @cvar PCG: The preconditioned conjugate gradient method (can only be applied
1217 for symmetric PDEs)
1218 @cvar CR: The conjugate residual method
1219 @cvar CGS: The conjugate gradient square method
1220 @cvar BICGSTAB: The stabilized BiConjugate Gradient method
1221 @cvar TFQMR: Transport Free Quasi Minimal Residual method
1222 @cvar MINRES: Minimum residual method
1223 @cvar SSOR: The symmetric overrelaxation method
1224 @cvar ILU0: The incomplete LU factorization preconditioner with no fill-in
1225 @cvar ILUT: The incomplete LU factorization preconditioner with fill-in
1226 @cvar JACOBI: The Jacobi preconditioner
1227 @cvar GMRES: The Gram-Schmidt minimum residual method
1228 @cvar PRES20: Special GMRES with restart after 20 steps and truncation after
1229 5 residuals
1230 @cvar LUMPING: Matrix lumping
1231 @cvar NO_REORDERING: No matrix reordering allowed
1232 @cvar MINIMUM_FILL_IN: Reorder matrix to reduce fill-in during factorization
1233 @cvar NESTED_DISSECTION: Reorder matrix to improve load balancing during
1234 factorization
1235 @cvar PASO: PASO solver package
1236 @cvar SCSL: SGI SCSL solver library
1237 @cvar MKL: Intel's MKL solver library
1238 @cvar UMFPACK: The UMFPACK library
1239 @cvar TRILINOS: The TRILINOS parallel solver class library from Sandia Natl
1240 Labs
1241 @cvar ITERATIVE: The default iterative solver
1242 @cvar AMG: Algebraic Multi Grid
1243 @cvar RILU: Recursive ILU
1244 @cvar GS: Gauss-Seidel solver
1245
1246 """
1247 DEFAULT= 0
1248 DIRECT= 1
1249 CHOLEVSKY= 2
1250 PCG= 3
1251 CR= 4
1252 CGS= 5
1253 BICGSTAB= 6
1254 SSOR= 7
1255 ILU0= 8
1256 ILUT= 9
1257 JACOBI= 10
1258 GMRES= 11
1259 PRES20= 12
1260 LUMPING= 13
1261 NO_REORDERING= 17
1262 MINIMUM_FILL_IN= 18
1263 NESTED_DISSECTION= 19
1264 SCSL= 14
1265 MKL= 15
1266 UMFPACK= 16
1267 ITERATIVE= 20
1268 PASO= 21
1269 AMG= 22
1270 RILU = 23
1271 TRILINOS = 24
1272 NONLINEAR_GMRES = 25
1273 TFQMR = 26
1274 MINRES = 27
1275 GS=28
1276
1277 SMALL_TOLERANCE=1.e-13
1278 PACKAGE_KEY="package"
1279 METHOD_KEY="method"
1280 SYMMETRY_KEY="symmetric"
1281 TOLERANCE_KEY="tolerance"
1282 PRECONDITIONER_KEY="preconditioner"
1283
1284
1285 def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):
1286 """
1287 Initializes a linear problem.
1288
1289 @param domain: domain of the PDE
1290 @type domain: L{Domain<escript.Domain>}
1291 @param numEquations: number of equations. If C{None} the number of
1292 equations is extracted from the coefficients.
1293 @param numSolutions: number of solution components. If C{None} the number
1294 of solution components is extracted from the
1295 coefficients.
1296 @param debug: if True debug information is printed
1297
1298 """
1299 super(LinearProblem, self).__init__()
1300
1301 self.__debug=debug
1302 self.__domain=domain
1303 self.__numEquations=numEquations
1304 self.__numSolutions=numSolutions
1305 self.__altered_coefficients=False
1306 self.__reduce_equation_order=False
1307 self.__reduce_solution_order=False
1308 self.__tolerance=1.e-8
1309 self.__solver_method=self.DEFAULT
1310 self.__solver_package=self.DEFAULT
1311 self.__preconditioner=self.DEFAULT
1312 self.__is_RHS_valid=False
1313 self.__is_operator_valid=False
1314 self.__sym=False
1315 self.__COEFFICIENTS={}
1316 # initialize things:
1317 self.resetAllCoefficients()
1318 self.__system_type=None
1319 self.updateSystemType()
1320 # ==========================================================================
1321 # general stuff:
1322 # ==========================================================================
1323 def __str__(self):
1324 """
1325 Returns a string representation of the PDE.
1326
1327 @return: a simple representation of the PDE
1328 @rtype: C{str}
1329 """
1330 return "<LinearProblem %d>"%id(self)
1331 # ==========================================================================
1332 # debug :
1333 # ==========================================================================
1334 def setDebugOn(self):
1335 """
1336 Switches debug output on.
1337 """
1338 self.__debug=not None
1339
1340 def setDebugOff(self):
1341 """
1342 Switches debug output off.
1343 """
1344 self.__debug=None
1345
1346 def setDebug(self, flag):
1347 """
1348 Switches debug output on if C{flag} is True otherwise it is switched off.
1349
1350 @param flag: desired debug status
1351 @type flag: C{bool}
1352 """
1353 if flag:
1354 self.setDebugOn()
1355 else:
1356 self.setDebugOff()
1357
1358 def trace(self,text):
1359 """
1360 Prints the text message if debug mode is switched on.
1361
1362 @param text: message to be printed
1363 @type text: C{string}
1364 """
1365 if self.__debug: print "%s: %s"%(str(self),text)
1366
1367 # ==========================================================================
1368 # some service functions:
1369 # ==========================================================================
1370 def introduceCoefficients(self,**coeff):
1371 """
1372 Introduces new coefficients into the problem.
1373
1374 Use::
1375
1376 p.introduceCoefficients(A=PDECoef(...), B=PDECoef(...))
1377
1378 to introduce the coefficients M{A} ans M{B}.
1379 """
1380 for name, type in coeff.items():
1381 if not isinstance(type,PDECoef):
1382 raise ValueError,"coefficient %s has no type."%name
1383 self.__COEFFICIENTS[name]=type
1384 self.__COEFFICIENTS[name].resetValue()
1385 self.trace("coefficient %s has been introduced."%name)
1386
1387 def getDomain(self):
1388 """
1389 Returns the domain of the PDE.
1390
1391 @return: the domain of the PDE
1392 @rtype: L{Domain<escript.Domain>}
1393 """
1394 return self.__domain
1395
1396 def getDim(self):
1397 """
1398 Returns the spatial dimension of the PDE.
1399
1400 @return: the spatial dimension of the PDE domain
1401 @rtype: C{int}
1402 """
1403 return self.getDomain().getDim()
1404
1405 def getNumEquations(self):
1406 """
1407 Returns the number of equations.
1408
1409 @return: the number of equations
1410 @rtype: C{int}
1411 @raise UndefinedPDEError: if the number of equations is not specified yet
1412 """
1413 if self.__numEquations==None:
1414 if self.__numSolutions==None:
1415 raise UndefinedPDEError,"Number of equations is undefined. Please specify argument numEquations."
1416 else:
1417 self.__numEquations=self.__numSolutions
1418 return self.__numEquations
1419
1420 def getNumSolutions(self):
1421 """
1422 Returns the number of unknowns.
1423
1424 @return: the number of unknowns
1425 @rtype: C{int}
1426 @raise UndefinedPDEError: if the number of unknowns is not specified yet
1427 """
1428 if self.__numSolutions==None:
1429 if self.__numEquations==None:
1430 raise UndefinedPDEError,"Number of solution is undefined. Please specify argument numSolutions."
1431 else:
1432 self.__numSolutions=self.__numEquations
1433 return self.__numSolutions
1434
1435 def reduceEquationOrder(self):
1436 """
1437 Returns the status of order reduction for the equation.
1438
1439 @return: True if reduced interpolation order is used for the
1440 representation of the equation, False otherwise
1441 @rtype: L{bool}
1442 """
1443 return self.__reduce_equation_order
1444
1445 def reduceSolutionOrder(self):
1446 """
1447 Returns the status of order reduction for the solution.
1448
1449 @return: True if reduced interpolation order is used for the
1450 representation of the solution, False otherwise
1451 @rtype: L{bool}
1452 """
1453 return self.__reduce_solution_order
1454
1455 def getFunctionSpaceForEquation(self):
1456 """
1457 Returns the L{FunctionSpace<escript.FunctionSpace>} used to discretize
1458 the equation.
1459
1460 @return: representation space of equation
1461 @rtype: L{FunctionSpace<escript.FunctionSpace>}
1462 """
1463 if self.reduceEquationOrder():
1464 return escript.ReducedSolution(self.getDomain())
1465 else:
1466 return escript.Solution(self.getDomain())
1467
1468 def getFunctionSpaceForSolution(self):
1469 """
1470 Returns the L{FunctionSpace<escript.FunctionSpace>} used to represent
1471 the solution.
1472
1473 @return: representation space of solution
1474 @rtype: L{FunctionSpace<escript.FunctionSpace>}
1475 """
1476 if self.reduceSolutionOrder():
1477 return escript.ReducedSolution(self.getDomain())
1478 else:
1479 return escript.Solution(self.getDomain())
1480
1481 # ==========================================================================
1482 # solver settings:
1483 # ==========================================================================
1484 def setSolverMethod(self,solver=None,preconditioner=None):
1485 """
1486 Sets a new solver method and/or preconditioner.
1487
1488 @param solver: new solver method to use
1489 @type solver: one of L{DEFAULT}, L{ITERATIVE} L{DIRECT}, L{CHOLEVSKY},
1490 L{PCG}, L{CR}, L{CGS}, L{BICGSTAB}, L{SSOR}, L{GMRES},
1491 L{TFQMR}, L{MINRES}, L{PRES20}, L{LUMPING}, L{AMG}
1492 @param preconditioner: new preconditioner to use
1493 @type preconditioner: one of L{DEFAULT}, L{JACOBI} L{ILU0}, L{ILUT},
1494 L{SSOR}, L{RILU}, L{GS}
1495
1496 @note: the solver method chosen may not be available by the selected
1497 package.
1498 """
1499 if solver==None: solver=self.__solver_method
1500 if preconditioner==None: preconditioner=self.__preconditioner
1501 if solver==None: solver=self.DEFAULT
1502 if preconditioner==None: preconditioner=self.DEFAULT
1503 if not (solver,preconditioner)==self.getSolverMethod():
1504 self.__solver_method=solver
1505 self.__preconditioner=preconditioner
1506 self.updateSystemType()
1507 self.trace("New solver is %s"%self.getSolverMethodName())
1508
1509 def getSolverMethodName(self):
1510 """
1511 Returns the name of the solver currently used.
1512
1513 @return: the name of the solver currently used
1514 @rtype: C{string}
1515 """
1516
1517 m=self.getSolverMethod()
1518 p=self.getSolverPackage()
1519 method=""
1520 if m[0]==self.DEFAULT: method="DEFAULT"
1521 elif m[0]==self.DIRECT: method= "DIRECT"
1522 elif m[0]==self.ITERATIVE: method= "ITERATIVE"
1523 elif m[0]==self.CHOLEVSKY: method= "CHOLEVSKY"
1524 elif m[0]==self.PCG: method= "PCG"
1525 elif m[0]==self.TFQMR: method= "TFQMR"
1526 elif m[0]==self.MINRES: method= "MINRES"
1527 elif m[0]==self.CR: method= "CR"
1528 elif m[0]==self.CGS: method= "CGS"
1529 elif m[0]==self.BICGSTAB: method= "BICGSTAB"
1530 elif m[0]==self.SSOR: method= "SSOR"
1531 elif m[0]==self.GMRES: method= "GMRES"
1532 elif m[0]==self.PRES20: method= "PRES20"
1533 elif m[0]==self.LUMPING: method= "LUMPING"
1534 elif m[0]==self.AMG: method= "AMG"
1535 if m[1]==self.DEFAULT: method+="+DEFAULT"
1536 elif m[1]==self.JACOBI: method+= "+JACOBI"
1537 elif m[1]==self.ILU0: method+= "+ILU0"
1538 elif m[1]==self.ILUT: method+= "+ILUT"
1539 elif m[1]==self.SSOR: method+= "+SSOR"
1540 elif m[1]==self.AMG: method+= "+AMG"
1541 elif m[1]==self.RILU: method+= "+RILU"
1542 elif m[1]==self.GS: method+= "+GS"
1543 if p==self.DEFAULT: package="DEFAULT"
1544 elif p==self.PASO: package= "PASO"
1545 elif p==self.MKL: package= "MKL"
1546 elif p==self.SCSL: package= "SCSL"
1547 elif p==self.UMFPACK: package= "UMFPACK"
1548 elif p==self.TRILINOS: package= "TRILINOS"
1549 else : method="unknown"
1550 return "%s solver of %s package"%(method,package)
1551
1552 def getSolverMethod(self):
1553 """
1554 Returns the solver method and preconditioner used.
1555
1556 @return: the solver method currently used.
1557 @rtype: C{tuple} of C{int}
1558 """
1559 return self.__solver_method,self.__preconditioner
1560
1561 def setSolverPackage(self,package=None):
1562 """
1563 Sets a new solver package.
1564
1565 @param package: the new solver package
1566 @type package: one of L{DEFAULT}, L{PASO} L{SCSL}, L{MKL}, L{UMFPACK},
1567 L{TRILINOS}
1568 """
1569 if package==None: package=self.DEFAULT
1570 if not package==self.getSolverPackage():
1571 self.__solver_package=package
1572 self.updateSystemType()
1573 self.trace("New solver is %s"%self.getSolverMethodName())
1574
1575 def getSolverPackage(self):
1576 """
1577 Returns the package of the solver.
1578
1579 @return: the solver package currently being used
1580 @rtype: C{int}
1581 """
1582 return self.__solver_package
1583
1584 def isUsingLumping(self):
1585 """
1586 Checks if matrix lumping is the current solver method.
1587
1588 @return: True if the current solver method is lumping
1589 @rtype: C{bool}
1590 """
1591 return self.getSolverMethod()[0]==self.LUMPING
1592
1593 def setTolerance(self,tol=1.e-8):
1594 """
1595 Resets the tolerance for the solver method to C{tol}.
1596
1597 @param tol: new tolerance for the solver. If C{tol} is lower than the
1598 current tolerence the system will be resolved.
1599 @type tol: positive C{float}
1600 @raise ValueError: if tolerance is not positive
1601 """
1602 if not tol>0:
1603 raise ValueError,"Tolerance has to be positive"
1604 if tol<self.getTolerance(): self.invalidateSolution()
1605 self.trace("New tolerance %e"%tol)
1606 self.__tolerance=tol
1607 return
1608
1609 def getTolerance(self):
1610 """
1611 Returns the tolerance currently set for the solution.
1612
1613 @return: tolerance currently used
1614 @rtype: C{float}
1615 """
1616 return self.__tolerance
1617
1618 # ==========================================================================
1619 # symmetry flag:
1620 # ==========================================================================
1621 def isSymmetric(self):
1622 """
1623 Checks if symmetry is indicated.
1624
1625 @return: True if a symmetric PDE is indicated, False otherwise
1626 @rtype: C{bool}
1627 """
1628 return self.__sym
1629
1630 def setSymmetryOn(self):
1631 """
1632 Sets the symmetry flag.
1633 """
1634 if not self.isSymmetric():
1635 self.trace("PDE is set to be symmetric")
1636 self.__sym=True
1637 self.updateSystemType()
1638
1639 def setSymmetryOff(self):
1640 """
1641 Clears the symmetry flag.
1642 """
1643 if self.isSymmetric():
1644 self.trace("PDE is set to be nonsymmetric")
1645 self.__sym=False
1646 self.updateSystemType()
1647
1648 def setSymmetryTo(self,flag=False):
1649 """
1650 Sets the symmetry flag to C{flag}.
1651
1652 @param flag: If True, the symmetry flag is set otherwise reset.
1653 @type flag: C{bool}
1654 """
1655 if flag:
1656 self.setSymmetryOn()
1657 else:
1658 self.setSymmetryOff()
1659
1660 # ==========================================================================
1661 # function space handling for the equation as well as the solution
1662 # ==========================================================================
1663 def setReducedOrderOn(self):
1664 """
1665 Switches reduced order on for solution and equation representation.
1666
1667 @raise RuntimeError: if order reduction is altered after a coefficient has
1668 been set
1669 """
1670 self.setReducedOrderForSolutionOn()
1671 self.setReducedOrderForEquationOn()
1672
1673 def setReducedOrderOff(self):
1674 """
1675 Switches reduced order off for solution and equation representation
1676
1677 @raise RuntimeError: if order reduction is altered after a coefficient has
1678 been set
1679 """
1680 self.setReducedOrderForSolutionOff()
1681 self.setReducedOrderForEquationOff()
1682
1683 def setReducedOrderTo(self,flag=False):
1684 """
1685 Sets order reduction state for both solution and equation representation
1686 according to flag.
1687
1688 @param flag: if True, the order reduction is switched on for both solution
1689 and equation representation, otherwise or if flag is not
1690 present order reduction is switched off
1691 @type flag: C{bool}
1692 @raise RuntimeError: if order reduction is altered after a coefficient has
1693 been set
1694 """
1695 self.setReducedOrderForSolutionTo(flag)
1696 self.setReducedOrderForEquationTo(flag)
1697
1698
1699 def setReducedOrderForSolutionOn(self):
1700 """
1701 Switches reduced order on for solution representation.
1702
1703 @raise RuntimeError: if order reduction is altered after a coefficient has
1704 been set
1705 """
1706 if not self.__reduce_solution_order:
1707 if self.__altered_coefficients:
1708 raise RuntimeError,"order cannot be altered after coefficients have been defined."
1709 self.trace("Reduced order is used for solution representation.")
1710 self.__reduce_solution_order=True
1711 self.initializeSystem()
1712
1713 def setReducedOrderForSolutionOff(self):
1714 """
1715 Switches reduced order off for solution representation
1716
1717 @raise RuntimeError: if order reduction is altered after a coefficient has
1718 been set.
1719 """
1720 if self.__reduce_solution_order:
1721 if self.__altered_coefficients:
1722 raise RuntimeError,"order cannot be altered after coefficients have been defined."
1723 self.trace("Full order is used to interpolate solution.")
1724 self.__reduce_solution_order=False
1725 self.initializeSystem()
1726
1727 def setReducedOrderForSolutionTo(self,flag=False):
1728 """
1729 Sets order reduction state for solution representation according to flag.
1730
1731 @param flag: if flag is True, the order reduction is switched on for
1732 solution representation, otherwise or if flag is not present
1733 order reduction is switched off
1734 @type flag: C{bool}
1735 @raise RuntimeError: if order reduction is altered after a coefficient has
1736 been set
1737 """
1738 if flag:
1739 self.setReducedOrderForSolutionOn()
1740 else:
1741 self.setReducedOrderForSolutionOff()
1742
1743 def setReducedOrderForEquationOn(self):
1744 """
1745 Switches reduced order on for equation representation.
1746
1747 @raise RuntimeError: if order reduction is altered after a coefficient has
1748 been set
1749 """
1750 if not self.__reduce_equation_order:
1751 if self.__altered_coefficients:
1752 raise RuntimeError,"order cannot be altered after coefficients have been defined."
1753 self.trace("Reduced order is used for test functions.")
1754 self.__reduce_equation_order=True
1755 self.initializeSystem()
1756
1757 def setReducedOrderForEquationOff(self):
1758 """
1759 Switches reduced order off for equation representation.
1760
1761 @raise RuntimeError: if order reduction is altered after a coefficient has
1762 been set
1763 """
1764 if self.__reduce_equation_order:
1765 if self.__altered_coefficients:
1766 raise RuntimeError,"order cannot be altered after coefficients have been defined."
1767 self.trace("Full order is used for test functions.")
1768 self.__reduce_equation_order=False
1769 self.initializeSystem()
1770
1771 def setReducedOrderForEquationTo(self,flag=False):
1772 """
1773 Sets order reduction state for equation representation according to flag.
1774
1775 @param flag: if flag is True, the order reduction is switched on for
1776 equation representation, otherwise or if flag is not present
1777 order reduction is switched off
1778 @type flag: C{bool}
1779 @raise RuntimeError: if order reduction is altered after a coefficient has
1780 been set
1781 """
1782 if flag:
1783 self.setReducedOrderForEquationOn()
1784 else:
1785 self.setReducedOrderForEquationOff()
1786
1787 def updateSystemType(self):
1788 """
1789 Reassesses the system type and, if a new matrix is needed, resets the
1790 system.
1791 """
1792 new_system_type=self.getRequiredSystemType()
1793 if not new_system_type==self.__system_type:
1794 self.trace("Matrix type is now %d."%new_system_type)
1795 self.__system_type=new_system_type
1796 self.initializeSystem()
1797
1798 def getSystemType(self):
1799 """
1800 Returns the current system type.
1801 """
1802 return self.__system_type
1803
1804 def checkSymmetricTensor(self,name,verbose=True):
1805 """
1806 Tests a coefficient for symmetry.
1807
1808 @param name: name of the coefficient
1809 @type name: C{str}
1810 @param verbose: if set to True or not present a report on coefficients
1811 which break the symmetry is printed.
1812 @type verbose: C{bool}
1813 @return: True if coefficient C{name} is symmetric
1814 @rtype: C{bool}
1815 """
1816 A=self.getCoefficient(name)
1817 verbose=verbose or self.__debug
1818 out=True
1819 if not A.isEmpty():
1820 tol=util.Lsup(A)*self.SMALL_TOLERANCE
1821 s=A.getShape()
1822 if A.getRank() == 4:
1823 if s[0]==s[2] and s[1] == s[3]:
1824 for i in range(s[0]):
1825 for j in range(s[1]):
1826 for k in range(s[2]):
1827 for l in range(s[3]):
1828 if util.Lsup(A[i,j,k,l]-A[k,l,i,j])>tol:
1829 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)
1830 out=False
1831 else:
1832 if verbose: print "non-symmetric problem because of inappropriate shape %s of coefficient %s."%(s,name)
1833 out=False
1834 elif A.getRank() == 2:
1835 if s[0]==s[1]:
1836 for j in range(s[0]):
1837 for l in range(s[1]):
1838 if util.Lsup(A[j,l]-A[l,j])>tol:
1839 if verbose: print "non-symmetric problem because %s[%d,%d]!=%s[%d,%d]"%(name,j,l,name,l,j)
1840 out=False
1841 else:
1842 if verbose: print "non-symmetric problem because of inappropriate shape %s of coefficient %s."%(s,name)
1843 out=False
1844 elif A.getRank() == 0:
1845 pass
1846 else:
1847 raise ValueError,"Cannot check rank %s of %s."%(A.getRank(),name)
1848 return out
1849
1850 def checkReciprocalSymmetry(self,name0,name1,verbose=True):
1851 """
1852 Tests two coefficients for reciprocal symmetry.
1853
1854 @param name0: name of the first coefficient
1855 @type name0: C{str}
1856 @param name1: name of the second coefficient
1857 @type name1: C{str}
1858 @param verbose: if set to True or not present a report on coefficients
1859 which break the symmetry is printed
1860 @type verbose: C{bool}
1861 @return: True if coefficients C{name0} and C{name1} are reciprocally
1862 symmetric.
1863 @rtype: C{bool}
1864 """
1865 B=self.getCoefficient(name0)
1866 C=self.getCoefficient(name1)
1867 verbose=verbose or self.__debug
1868 out=True
1869 if B.isEmpty() and not C.isEmpty():
1870 if verbose: print "non-symmetric problem because %s is not present but %s is"%(name0,name1)
1871 out=False
1872 elif not B.isEmpty() and C.isEmpty():
1873 if verbose: print "non-symmetric problem because %s is not present but %s is"%(name0,name1)
1874 out=False
1875 elif not B.isEmpty() and not C.isEmpty():
1876 sB=B.getShape()
1877 sC=C.getShape()
1878 tol=(util.Lsup(B)+util.Lsup(C))*self.SMALL_TOLERANCE/2.
1879 if len(sB) != len(sC):
1880 if verbose: print "non-symmetric problem because ranks of %s (=%s) and %s (=%s) are different."%(name0,len(sB),name1,len(sC))
1881 out=False
1882 else:
1883 if len(sB)==0:
1884 if util.Lsup(B-C)>tol:
1885 if verbose: print "non-symmetric problem because %s!=%s"%(name0,name1)
1886 out=False
1887 elif len(sB)==1:
1888 if sB[0]==sC[0]:
1889 for j in range(sB[0]):
1890 if util.Lsup(B[j]-C[j])>tol:
1891 if verbose: print "non-symmetric PDE because %s[%d]!=%s[%d]"%(name0,j,name1,j)
1892 out=False
1893 else:
1894 if verbose: print "non-symmetric problem because of inappropriate shapes %s and %s of coefficients %s and %s, respectively."%(sB,sC,name0,name1)
1895 elif len(sB)==3:
1896 if sB[0]==sC[1] and sB[1]==sC[2] and sB[2]==sC[0]:
1897 for i in range(sB[0]):
1898 for j in range(sB[1]):
1899 for k in range(sB[2]):
1900 if util.Lsup(B[i,j,k]-C[k,i,j])>tol:
1901 if verbose: print "non-symmetric problem because %s[%d,%d,%d]!=%s[%d,%d,%d]"%(name0,i,j,k,name1,k,i,j)
1902 out=False
1903 else:
1904 if verbose: print "non-symmetric problem because of inappropriate shapes %s and %s of coefficients %s and %s, respectively."%(sB,sC,name0,name1)
1905 else:
1906 raise ValueError,"Cannot check rank %s of %s and %s."%(len(sB),name0,name1)
1907 return out
1908
1909 def getCoefficient(self,name):
1910 """
1911 Returns the value of the coefficient C{name}.
1912
1913 @param name: name of the coefficient requested
1914 @type name: C{string}
1915 @return: the value of the coefficient
1916 @rtype: L{Data<escript.Data>}
1917 @raise IllegalCoefficient: if C{name} is not a coefficient of the PDE
1918 """
1919 if self.hasCoefficient(name):
1920 return self.__COEFFICIENTS[name].getValue()
1921 else:
1922 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1923
1924 def hasCoefficient(self,name):
1925 """
1926 Returns True if C{name} is the name of a coefficient.
1927
1928 @param name: name of the coefficient enquired
1929 @type name: C{string}
1930 @return: True if C{name} is the name of a coefficient of the general PDE,
1931 False otherwise
1932 @rtype: C{bool}
1933 """
1934 return self.__COEFFICIENTS.has_key(name)
1935
1936 def createCoefficient(self, name):
1937 """
1938 Creates a L{Data<escript.Data>} object corresponding to coefficient
1939 C{name}.
1940
1941 @return: the coefficient C{name} initialized to 0
1942 @rtype: L{Data<escript.Data>}
1943 @raise IllegalCoefficient: if C{name} is not a coefficient of the PDE
1944 """
1945 if self.hasCoefficient(name):
1946 return escript.Data(0.,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name))
1947 else:
1948 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1949
1950 def getFunctionSpaceForCoefficient(self,name):
1951 """
1952 Returns the L{FunctionSpace<escript.FunctionSpace>} to be used for
1953 coefficient C{name}.
1954
1955 @param name: name of the coefficient enquired
1956 @type name: C{string}
1957 @return: the function space to be used for coefficient C{name}
1958 @rtype: L{FunctionSpace<escript.FunctionSpace>}
1959 @raise IllegalCoefficient: if C{name} is not a coefficient of the PDE
1960 """
1961 if self.hasCoefficient(name):
1962 return self.__COEFFICIENTS[name].getFunctionSpace(self.getDomain())
1963 else:
1964 raise ValueError,"unknown coefficient %s requested"%name
1965
1966 def getShapeOfCoefficient(self,name):
1967 """
1968 Returns the shape of the coefficient C{name}.
1969
1970 @param name: name of the coefficient enquired
1971 @type name: C{string}
1972 @return: the shape of the coefficient C{name}
1973 @rtype: C{tuple} of C{int}
1974 @raise IllegalCoefficient: if C{name} is not a coefficient of the PDE
1975 """
1976 if self.hasCoefficient(name):
1977 return self.__COEFFICIENTS[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())
1978 else:
1979 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1980
1981 def resetAllCoefficients(self):
1982 """
1983 Resets all coefficients to their default values.
1984 """
1985 for i in self.__COEFFICIENTS.iterkeys():
1986 self.__COEFFICIENTS[i].resetValue()
1987
1988 def alteredCoefficient(self,name):
1989 """
1990 Announces that coefficient C{name} has been changed.
1991
1992 @param name: name of the coefficient affected
1993 @type name: C{string}
1994 @raise IllegalCoefficient: if C{name} is not a coefficient of the PDE
1995 @note: if C{name} is q or r, the method will not trigger a rebuild of the
1996 system as constraints are applied to the solved system.
1997 """
1998 if self.hasCoefficient(name):
1999 self.trace("Coefficient %s has been altered."%name)
2000 if not ((name=="q" or name=="r") and self.isUsingLumping()):
2001 if self.__COEFFICIENTS[name].isAlteringOperator(): self.invalidateOperator()
2002 if self.__COEFFICIENTS[name].isAlteringRightHandSide(): self.invalidateRightHandSide()
2003 else:
2004 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
2005
2006 def validSolution(self):
2007 """
2008 Marks the solution as valid.
2009 """
2010 self.__is_solution_valid=True
2011
2012 def invalidateSolution(self):
2013 """
2014 Indicates the PDE has to be resolved if the solution is requested.
2015 """
2016 self.trace("System will be resolved.")
2017 self.__is_solution_valid=False
2018
2019 def isSolutionValid(self):
2020 """
2021 Returns True if the solution is still valid.
2022 """
2023 return self.__is_solution_valid
2024
2025 def validOperator(self):
2026 """
2027 Marks the operator as valid.
2028 """
2029 self.__is_operator_valid=True
2030
2031 def invalidateOperator(self):
2032 """
2033 Indicates the operator has to be rebuilt next time it is used.
2034 """
2035 self.trace("Operator will be rebuilt.")
2036 self.invalidateSolution()
2037 self.__is_operator_valid=False
2038
2039 def isOperatorValid(self):
2040 """
2041 Returns True if the operator is still valid.
2042 """
2043 return self.__is_operator_valid
2044
2045 def validRightHandSide(self):
2046 """
2047 Marks the right hand side as valid.
2048 """
2049 self.__is_RHS_valid=True
2050
2051 def invalidateRightHandSide(self):
2052 """
2053 Indicates the right hand side has to be rebuilt next time it is used.
2054 """
2055 if self.isRightHandSideValid(): self.trace("Right hand side has to be rebuilt.")
2056 self.invalidateSolution()
2057 self.__is_RHS_valid=False
2058
2059 def isRightHandSideValid(self):
2060 """
2061 Returns True if the operator is still valid.
2062 """
2063 return self.__is_RHS_valid
2064
2065 def invalidateSystem(self):
2066 """
2067 Announces that everything has to be rebuilt.
2068 """
2069 if self.isRightHandSideValid(): self.trace("System has to be rebuilt.")
2070 self.invalidateSolution()
2071 self.invalidateOperator()
2072 self.invalidateRightHandSide()
2073
2074 def isSystemValid(self):
2075 """
2076 Returns True if the system (including solution) is still vaild.
2077 """
2078 return self.isSolutionValid() and self.isOperatorValid() and self.isRightHandSideValid()
2079
2080 def initializeSystem(self):
2081 """
2082 Resets the system clearing the operator, right hand side and solution.
2083 """
2084 self.trace("New System has been created.")
2085 self.__operator=escript.Operator()
2086 self.__righthandside=escript.Data()
2087 self.__solution=escript.Data()
2088 self.invalidateSystem()
2089
2090 def getOperator(self):
2091 """
2092 Returns the operator of the linear problem.
2093
2094 @return: the operator of the problem
2095 """
2096 return self.getSystem()[0]
2097
2098 def getRightHandSide(self):
2099 """
2100 Returns the right hand side of the linear problem.
2101
2102 @return: the right hand side of the problem
2103 @rtype: L{Data<escript.Data>}
2104 """
2105 return self.getSystem()[1]
2106
2107 def createRightHandSide(self):
2108 """
2109 Returns an instance of a new right hand side.
2110 """
2111 self.trace("New right hand side is allocated.")
2112 if self.getNumEquations()>1:
2113 return escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True)
2114 else:
2115 return escript.Data(0.,(),self.getFunctionSpaceForEquation(),True)
2116
2117 def createSolution(self):
2118 """
2119 Returns an instance of a new solution.
2120 """
2121 self.trace("New solution is allocated.")
2122 if self.getNumSolutions()>1:
2123 return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)
2124 else:
2125 return escript.Data(0.,(),self.getFunctionSpaceForSolution(),True)
2126
2127 def resetSolution(self):
2128 """
2129 Sets the solution to zero.
2130 """
2131 if self.__solution.isEmpty():
2132 self.__solution=self.createSolution()
2133 else:
2134 self.__solution.setToZero()
2135 self.trace("Solution is reset to zero.")
2136
2137 def setSolution(self,u):
2138 """
2139 Sets the solution assuming that makes the system valid.
2140 """
2141 self.__solution=u
2142 self.validSolution()
2143
2144 def getCurrentSolution(self):
2145 """
2146 Returns the solution in its current state.
2147 """
2148 if self.__solution.isEmpty(): self.__solution=self.createSolution()
2149 return self.__solution
2150
2151 def resetRightHandSide(self):
2152 """
2153 Sets the right hand side to zero.
2154 """
2155 if self.__righthandside.isEmpty():
2156 self.__righthandside=self.createRightHandSide()
2157 else:
2158 self.__righthandside.setToZero()
2159 self.trace("Right hand side is reset to zero.")
2160
2161 def getCurrentRightHandSide(self):
2162 """
2163 Returns the right hand side in its current state.
2164 """
2165 if self.__righthandside.isEmpty(): self.__righthandside=self.createRightHandSide()
2166 return self.__righthandside
2167
2168 def resetOperator(self):
2169 """
2170 Makes sure that the operator is instantiated and returns it initialized
2171 with zeros.
2172 """
2173 if self.__operator.isEmpty():
2174 if self.isUsingLumping():
2175 self.__operator=self.createSolution()
2176 else:
2177 self.__operator=self.createOperator()
2178 else:
2179 if self.isUsingLumping():
2180 self.__operator.setToZero()
2181 else:
2182 self.__operator.resetValues()
2183 self.trace("Operator reset to zero")
2184
2185 def getCurrentOperator(self):
2186 """
2187 Returns the operator in its current state.
2188 """
2189 if self.__operator.isEmpty(): self.resetOperator()
2190 return self.__operator
2191
2192 def setValue(self,**coefficients):
2193 """
2194 Sets new values to coefficients.
2195
2196 @raise IllegalCoefficient: if an unknown coefficient keyword is used
2197 """
2198 # check if the coefficients are legal:
2199 for i in coefficients.iterkeys():
2200 if not self.hasCoefficient(i):
2201 raise IllegalCoefficient,"Attempt to set unknown coefficient %s"%i
2202 # if the number of unknowns or equations is still unknown we try to estimate them:
2203 if self.__numEquations==None or self.__numSolutions==None:
2204 for i,d in coefficients.iteritems():
2205 if hasattr(d,"shape"):
2206 s=d.shape
2207 elif hasattr(d,"getShape"):
2208 s=d.getShape()
2209 else:
2210 s=numpy.array(d).shape
2211 if s!=None:
2212 # get number of equations and number of unknowns:
2213 res=self.__COEFFICIENTS[i].estimateNumEquationsAndNumSolutions(self.getDomain(),s)
2214 if res==None:
2215 raise IllegalCoefficientValue,"Illegal shape %s of coefficient %s"%(s,i)
2216 else:
2217 if self.__numEquations==None: self.__numEquations=res[0]
2218 if self.__numSolutions==None: self.__numSolutions=res[1]
2219 if self.__numEquations==None: raise UndefinedPDEError,"unidentified number of equations"
2220 if self.__numSolutions==None: raise UndefinedPDEError,"unidentified number of solutions"
2221 # now we check the shape of the coefficient if numEquations and numSolutions are set:
2222 for i,d in coefficients.iteritems():
2223 try:
2224 self.__COEFFICIENTS[i].setValue(self.getDomain(),
2225 self.getNumEquations(),self.getNumSolutions(),
2226 self.reduceEquationOrder(),self.reduceSolutionOrder(),d)
2227 self.alteredCoefficient(i)
2228 except IllegalCoefficientFunctionSpace,m:
2229 # if the function space is wrong then we try the reduced version:
2230 i_red=i+"_reduced"
2231 if (not i_red in coefficients.keys()) and i_red in self.__COEFFICIENTS.keys():
2232 try:
2233 self.__COEFFICIENTS[i_red].setValue(self.getDomain(),
2234 self.getNumEquations(),self.getNumSolutions(),
2235 self.reduceEquationOrder(),self.reduceSolutionOrder(),d)
2236 self.alteredCoefficient(i_red)
2237 except IllegalCoefficientValue,m:
2238 raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
2239 except IllegalCoefficientFunctionSpace,m:
2240 raise IllegalCoefficientFunctionSpace("Coefficient %s:%s"%(i,m))
2241 else:
2242 raise IllegalCoefficientFunctionSpace("Coefficient %s:%s"%(i,m))
2243 except IllegalCoefficientValue,m:
2244 raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
2245 self.__altered_coefficients=True
2246
2247 # ==========================================================================
2248 # methods that are typically overwritten when implementing a particular
2249 # linear problem
2250 # ==========================================================================
2251 def getRequiredSystemType(self):
2252 """
2253 Returns the system type which needs to be used by the current set up.
2254
2255 @note: Typically this method is overwritten when implementing a
2256 particular linear problem.
2257 """
2258 return 0
2259
2260 def createOperator(self):
2261 """
2262 Returns an instance of a new operator.
2263
2264 @note: This method is overwritten when implementing a particular
2265 linear problem.
2266 """
2267 return escript.Operator()
2268
2269 def checkSymmetry(self,verbose=True):
2270 """
2271 Tests the PDE for symmetry.
2272
2273 @param verbose: if set to True or not present a report on coefficients
2274 which break the symmetry is printed
2275 @type verbose: C{bool}
2276 @return: True if the problem is symmetric
2277 @rtype: C{bool}
2278 @note: Typically this method is overwritten when implementing a
2279 particular linear problem.
2280 """
2281 out=True
2282 return out
2283
2284 def getSolution(self,**options):
2285 """
2286 Returns the solution of the problem.
2287
2288 @return: the solution
2289 @rtype: L{Data<escript.Data>}
2290
2291 @note: This method is overwritten when implementing a particular
2292 linear problem.
2293 """
2294 return self.getCurrentSolution()
2295
2296 def getSystem(self):
2297 """
2298 Returns the operator and right hand side of the PDE.
2299
2300 @return: the discrete version of the PDE
2301 @rtype: C{tuple} of L{Operator,<escript.Operator>} and L{Data<escript.Data>}.
2302
2303 @note: This method is overwritten when implementing a particular
2304 linear problem.
2305 """
2306 return (self.getCurrentOperator(), self.getCurrentRightHandSide())
2307 #=====================
2308
2309 class LinearPDE(LinearProblem):
2310 """
2311 This class is used to define a general linear, steady, second order PDE
2312 for an unknown function M{u} on a given domain defined through a
2313 L{Domain<escript.Domain>} object.
2314
2315 For a single PDE having a solution with a single component the linear PDE
2316 is defined in the following form:
2317
2318 M{-(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)}
2319
2320 where M{grad(F)} denotes the spatial derivative of M{F}. Einstein's
2321 summation convention, ie. summation over indexes appearing twice in a term
2322 of a sum performed, is used.
2323 The coefficients M{A}, M{B}, M{C}, M{D}, M{X} and M{Y} have to be specified
2324 through L{Data<escript.Data>} objects in L{Function<escript.Function>} and
2325 the coefficients M{A_reduced}, M{B_reduced}, M{C_reduced}, M{D_reduced},
2326 M{X_reduced} and M{Y_reduced} have to be specified through
2327 L{Data<escript.Data>} objects in L{ReducedFunction<escript.ReducedFunction>}.
2328 It is also allowed to use objects that can be converted into such
2329 L{Data<escript.Data>} objects. M{A} and M{A_reduced} are rank two, M{B},
2330 M{C}, M{X}, M{B_reduced}, M{C_reduced} and M{X_reduced} are rank one and
2331 M{D}, M{D_reduced}, M{Y} and M{Y_reduced} are scalar.
2332
2333 The following natural boundary conditions are considered:
2334
2335 M{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}
2336
2337 where M{n} is the outer normal field. Notice that the coefficients M{A},
2338 M{A_reduced}, M{B}, M{B_reduced}, M{X} and M{X_reduced} are defined in the
2339 PDE. The coefficients M{d} and M{y} are each a scalar in
2340 L{FunctionOnBoundary<escript.FunctionOnBoundary>} and the coefficients
2341 M{d_reduced} and M{y_reduced} are each a scalar in
2342 L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
2343
2344 Constraints for the solution prescribe the value of the solution at certain
2345 locations in the domain. They have the form
2346
2347 M{u=r} where M{q>0}
2348
2349 M{r} and M{q} are each scalar where M{q} is the characteristic function
2350 defining where the constraint is applied. The constraints override any
2351 other condition set by the PDE or the boundary condition.
2352
2353 The PDE is symmetrical if
2354
2355 M{A[i,j]=A[j,i]} and M{B[j]=C[j]} and M{A_reduced[i,j]=A_reduced[j,i]}
2356 and M{B_reduced[j]=C_reduced[j]}
2357
2358 For a system of PDEs and a solution with several components the PDE has the
2359 form
2360
2361 M{-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] }
2362
2363 M{A} and M{A_reduced} are of rank four, M{B}, M{B_reduced}, M{C} and
2364 M{C_reduced} are each of rank three, M{D}, M{D_reduced}, M{X_reduced} and
2365 M{X} are each of rank two and M{Y} and M{Y_reduced} are of rank one.
2366 The natural boundary conditions take the form:
2367
2368 M{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]}
2369
2370 The coefficient M{d} is of rank two and M{y} is of rank one both in
2371 L{FunctionOnBoundary<escript.FunctionOnBoundary>}. The coefficients
2372 M{d_reduced} is of rank two and M{y_reduced} is of rank one both in
2373 L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
2374
2375 Constraints take the form
2376
2377 M{u[i]=r[i]} where M{q[i]>0}
2378
2379 M{r} and M{q} are each rank one. Notice that at some locations not
2380 necessarily all components must have a constraint.
2381
2382 The system of PDEs is symmetrical if
2383
2384 - M{A[i,j,k,l]=A[k,l,i,j]}
2385 - M{A_reduced[i,j,k,l]=A_reduced[k,l,i,j]}
2386 - M{B[i,j,k]=C[k,i,j]}
2387 - M{B_reduced[i,j,k]=C_reduced[k,i,j]}
2388 - M{D[i,k]=D[i,k]}
2389 - M{D_reduced[i,k]=D_reduced[i,k]}
2390 - M{d[i,k]=d[k,i]}
2391 - M{d_reduced[i,k]=d_reduced[k,i]}
2392
2393 L{LinearPDE} also supports solution discontinuities over a contact region
2394 in the domain. To specify the conditions across the discontinuity we are
2395 using the generalised flux M{J} which, in the case of a system of PDEs
2396 and several components of the solution, is defined as
2397
2398 M{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]}
2399
2400 For the case of single solution component and single PDE M{J} is defined as
2401
2402 M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[j]+(B[i]+B_reduced[i])*u-X[i]-X_reduced[i]}
2403
2404 In the context of discontinuities M{n} denotes the normal on the
2405 discontinuity pointing from side 0 towards side 1 calculated from
2406 L{getNormal<escript.FunctionSpace.getNormal>} of L{FunctionOnContactZero<escript.FunctionOnContactZero>}.
2407 For a system of PDEs the contact condition takes the form
2408
2409 M{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]}
2410
2411 where M{J0} and M{J1} are the fluxes on side 0 and side 1 of the
2412 discontinuity, respectively. M{jump(u)}, which is the difference of the
2413 solution at side 1 and at side 0, denotes the jump of M{u} across
2414 discontinuity along the normal calculated by L{jump<util.jump>}.
2415 The coefficient M{d_contact} is of rank two and M{y_contact} is of rank one
2416 both in L{FunctionOnContactZero<escript.FunctionOnContactZero>} or
2417 L{FunctionOnContactOne<escript.FunctionOnContactOne>}.
2418 The coefficient M{d_contact_reduced} is of rank two and M{y_contact_reduced}
2419 is of rank one both in L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>}
2420 or L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}.
2421 In case of a single PDE and a single component solution the contact
2422 condition takes the form
2423
2424 M{n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)}
2425
2426 In this case the coefficient M{d_contact} and M{y_contact} are each scalar
2427 both in L{FunctionOnContactZero<escript.FunctionOnContactZero>} or
2428 L{FunctionOnContactOne<escript.FunctionOnContactOne>} and the coefficient
2429 M{d_contact_reduced} and M{y_contact_reduced} are each scalar both in
2430 L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or
2431 L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}.
2432
2433 Typical usage::
2434
2435 p = LinearPDE(dom)
2436 p.setValue(A=kronecker(dom), D=1, Y=0.5)
2437 u = p.getSolution()
2438
2439 """
2440
2441 def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):
2442 """
2443 Initializes a new linear PDE.
2444
2445 @param domain: domain of the PDE
2446 @type domain: L{Domain<escript.Domain>}
2447 @param numEquations: number of equations. If C{None} the number of
2448 equations is extracted from the PDE coefficients.
2449 @param numSolutions: number of solution components. If C{None} the number
2450 of solution components is extracted from the PDE
2451 coefficients.
2452 @param debug: if True debug information is printed
2453
2454 """
2455 super(LinearPDE, self).__init__(domain,numEquations,numSolutions,debug)
2456 #
2457 # the coefficients of the PDE:
2458 #
2459 self.introduceCoefficients(
2460 A=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2461 B=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2462 C=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2463 D=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2464 X=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
2465 Y=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2466 d=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2467 y=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2468 d_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2469 y_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2470 A_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2471 B_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2472 C_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2473 D_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2474 X_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
2475 Y_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2476 d_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2477 y_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2478 d_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2479 y_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2480 r=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.RIGHTHANDSIDE),
2481 q=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.BOTH) )
2482
2483 def __str__(self):
2484 """
2485 Returns the string representation of the PDE.
2486
2487 @return: a simple representation of the PDE
2488 @rtype: C{str}
2489 """
2490 return "<LinearPDE %d>"%id(self)
2491
2492 def getRequiredSystemType(self):
2493 """
2494 Returns the system type which needs to be used by the current set up.
2495 """
2496 return self.getDomain().getSystemMatrixTypeId(self.getSolverMethod()[0],self.getSolverMethod()[1],self.getSolverPackage(),self.isSymmetric())
2497
2498 def checkSymmetry(self,verbose=True):
2499 """
2500 Tests the PDE for symmetry.
2501
2502 @param verbose: if set to True or not present a report on coefficients
2503 which break the symmetry is printed.
2504 @type verbose: C{bool}
2505 @return: True if the PDE is symmetric
2506 @rtype: L{bool}
2507 @note: This is a very expensive operation. It should be used for
2508 degugging only! The symmetry flag is not altered.
2509 """
2510 out=True
2511 out=out and self.checkSymmetricTensor("A", verbose)
2512 out=out and self.checkSymmetricTensor("A_reduced", verbose)
2513 out=out and self.checkReciprocalSymmetry("B","C", verbose)
2514 out=out and self.checkReciprocalSymmetry("B_reduced","C_reduced", verbose)
2515 out=out and self.checkSymmetricTensor("D", verbose)
2516 out=out and self.checkSymmetricTensor("D_reduced", verbose)
2517 out=out and self.checkSymmetricTensor("d", verbose)
2518 out=out and self.checkSymmetricTensor("d_reduced", verbose)
2519 out=out and self.checkSymmetricTensor("d_contact", verbose)
2520 out=out and self.checkSymmetricTensor("d_contact_reduced", verbose)
2521 return out
2522
2523 def createOperator(self):
2524 """
2525 Returns an instance of a new operator.
2526 """
2527 self.trace("New operator is allocated.")
2528 return self.getDomain().newOperator( \
2529 self.getNumEquations(), \
2530 self.getFunctionSpaceForEquation(), \
2531 self.getNumSolutions(), \
2532 self.getFunctionSpaceForSolution(), \
2533 self.getSystemType())
2534
2535 def getSolution(self,**options):
2536 """
2537 Returns the solution of the PDE.
2538
2539 @return: the solution
2540 @rtype: L{Data<escript.Data>}
2541 @param options: solver options
2542 @keyword verbose: True to get some information during PDE solution
2543 @type verbose: C{bool}
2544 @keyword reordering: reordering scheme to be used during elimination.
2545 Allowed values are L{NO_REORDERING},
2546 L{MINIMUM_FILL_IN} and L{NESTED_DISSECTION}
2547 @keyword iter_max: maximum number of iteration steps allowed
2548 @keyword drop_tolerance: threshold for dropping in L{ILUT}
2549 @keyword drop_storage: maximum of allowed memory in L{ILUT}
2550 @keyword truncation: maximum number of residuals in L{GMRES}
2551 @keyword restart: restart cycle length in L{GMRES}
2552 """
2553 if not self.isSolutionValid():
2554 mat,f=self.getSystem()
2555 if self.isUsingLumping():
2556 self.setSolution(f*1/mat)
2557 else:
2558 options[self.TOLERANCE_KEY]=self.getTolerance()
2559 options[self.METHOD_KEY]=self.getSolverMethod()[0]
2560 options[self.PRECONDITIONER_KEY]=self.getSolverMethod()[1]
2561 options[self.PACKAGE_KEY]=self.getSolverPackage()
2562 options[self.SYMMETRY_KEY]=self.isSymmetric()
2563 self.trace("PDE is resolved.")
2564 self.trace("solver options: %s"%str(options))
2565 self.setSolution(mat.solve(f,options))
2566 return self.getCurrentSolution()
2567
2568 def getSystem(self):
2569 """
2570 Returns the operator and right hand side of the PDE.
2571
2572 @return: the discrete version of the PDE
2573 @rtype: C{tuple} of L{Operator,<escript.Operator>} and
2574 L{Data<escript.Data>}
2575 """
2576 if not self.isOperatorValid() or not self.isRightHandSideValid():
2577 if self.isUsingLumping():
2578 if not self.isOperatorValid():
2579 if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():
2580 raise TypeError,"Lumped matrix requires same order for equations and unknowns"
2581 if not self.getCoefficient("A").isEmpty():
2582 raise ValueError,"coefficient A in lumped matrix may not be present."
2583 if not self.getCoefficient("B").isEmpty():
2584 raise ValueError,"coefficient B in lumped matrix may not be present."
2585 if not self.getCoefficient("C").isEmpty():
2586 raise ValueError,"coefficient C in lumped matrix may not be present."
2587 if not self.getCoefficient("d_contact").isEmpty():
2588 raise ValueError,"coefficient d_contact in lumped matrix may not be present."
2589 if not self.getCoefficient("A_reduced").isEmpty():
2590 raise ValueError,"coefficient A_reduced in lumped matrix may not be present."
2591 if not self.getCoefficient("B_reduced").isEmpty():
2592 raise ValueError,"coefficient B_reduced in lumped matrix may not be present."
2593 if not self.getCoefficient("C_reduced").isEmpty():
2594 raise ValueError,"coefficient C_reduced in lumped matrix may not be present."
2595 if not self.getCoefficient("d_contact_reduced").isEmpty():
2596 raise ValueError,"coefficient d_contact_reduced in lumped matrix may not be present."
2597 D=self.getCoefficient("D")
2598 d=self.getCoefficient("d")
2599 D_reduced=self.getCoefficient("D_reduced")
2600 d_reduced=self.getCoefficient("d_reduced")
2601 if not D.isEmpty():
2602 if self.getNumSolutions()>1:
2603 D_times_e=util.matrix_mult(D,numpy.ones((self.getNumSolutions(),)))
2604 else:
2605 D_times_e=D
2606 else:
2607 D_times_e=escript.Data()
2608 if not d.isEmpty():
2609 if self.getNumSolutions()>1:
2610 d_times_e=util.matrix_mult(d,numpy.ones((self.getNumSolutions(),)))
2611 else:
2612 d_times_e=d
2613 else:
2614 d_times_e=escript.Data()
2615
2616 if not D_reduced.isEmpty():
2617 if self.getNumSolutions()>1:
2618 D_reduced_times_e=util.matrix_mult(D_reduced,numpy.ones((self.getNumSolutions(),)))
2619 else:
2620 D_reduced_times_e=D_reduced
2621 else:
2622 D_reduced_times_e=escript.Data()
2623 if not d_reduced.isEmpty():
2624 if self.getNumSolutions()>1:
2625 d_reduced_times_e=util.matrix_mult(d_reduced,numpy.ones((self.getNumSolutions(),)))
2626 else:
2627 d_reduced_times_e=d_reduced
2628 else:
2629 d_reduced_times_e=escript.Data()
2630
2631 self.resetOperator()
2632 operator=self.getCurrentOperator()
2633 if False and hasattr(self.getDomain(), "addPDEToLumpedSystem") :
2634 self.getDomain().addPDEToLumpedSystem(operator, D_times_e, d_times_e)
2635 self.getDomain().addPDEToLumpedSystem(operator, D_reduced_times_e, d_reduced_times_e)
2636 else:
2637 self.getDomain().addPDEToRHS(operator, \
2638 escript.Data(), \
2639 D_times_e, \
2640 d_times_e,\
2641 escript.Data())
2642 self.getDomain().addPDEToRHS(operator, \
2643 escript.Data(), \
2644 D_reduced_times_e, \
2645 d_reduced_times_e,\
2646 escript.Data())
2647 self.trace("New lumped operator has been built.")
2648 if not self.isRightHandSideValid():
2649 self.resetRightHandSide()
2650 righthandside=self.getCurrentRightHandSide()
2651 self.getDomain().addPDEToRHS(righthandside, \
2652 self.getCoefficient("X"), \
2653 self.getCoefficient("Y"),\
2654 self.getCoefficient("y"),\
2655 self.getCoefficient("y_contact"))
2656 self.getDomain().addPDEToRHS(righthandside, \
2657 self.getCoefficient("X_reduced"), \
2658 self.getCoefficient("Y_reduced"),\
2659 self.getCoefficient("y_reduced"),\
2660 self.getCoefficient("y_contact_reduced"))
2661 self.trace("New right hand side has been built.")
2662 self.validRightHandSide()
2663 self.insertConstraint(rhs_only=False)
2664 self.validOperator()
2665 else:
2666 if not self.isOperatorValid() and not self.isRightHandSideValid():
2667 self.resetRightHandSide()
2668 righthandside=self.getCurrentRightHandSide()
2669 self.resetOperator()
2670 operator=self.getCurrentOperator()
2671 self.getDomain().addPDEToSystem(operator,righthandside, \
2672 self.getCoefficient("A"), \
2673 self.getCoefficient("B"), \
2674 self.getCoefficient("C"), \
2675 self.getCoefficient("D"), \
2676 self.getCoefficient("X"), \
2677 self.getCoefficient("Y"), \
2678 self.getCoefficient("d"), \
2679 self.getCoefficient("y"), \
2680 self.getCoefficient("d_contact"), \
2681 self.getCoefficient("y_contact"))
2682 self.getDomain().addPDEToSystem(operator,righthandside, \
2683 self.getCoefficient("A_reduced"), \
2684 self.getCoefficient("B_reduced"), \
2685 self.getCoefficient("C_reduced"), \
2686 self.getCoefficient("D_reduced"), \
2687 self.getCoefficient("X_reduced"), \
2688 self.getCoefficient("Y_reduced"), \
2689 self.getCoefficient("d_reduced"), \
2690 self.getCoefficient("y_reduced"), \
2691 self.getCoefficient("d_contact_reduced"), \
2692 self.getCoefficient("y_contact_reduced"))
2693 self.insertConstraint(rhs_only=False)
2694 self.trace("New system has been built.")
2695 self.validOperator()
2696 self.validRightHandSide()
2697 elif not self.isRightHandSideValid():
2698 self.resetRightHandSide()
2699 righthandside=self.getCurrentRightHandSide()
2700 self.getDomain().addPDEToRHS(righthandside,
2701 self.getCoefficient("X"), \
2702 self.getCoefficient("Y"),\
2703 self.getCoefficient("y"),\
2704 self.getCoefficient("y_contact"))
2705 self.getDomain().addPDEToRHS(righthandside,
2706 self.getCoefficient("X_reduced"), \
2707 self.getCoefficient("Y_reduced"),\
2708 self.getCoefficient("y_reduced"),\
2709 self.getCoefficient("y_contact_reduced"))
2710 self.insertConstraint(rhs_only=True)
2711 self.trace("New right hand side has been built.")
2712 self.validRightHandSide()
2713 elif not self.isOperatorValid():
2714 self.resetOperator()
2715 operator=self.getCurrentOperator()
2716 self.getDomain().addPDEToSystem(operator,escript.Data(), \
2717 self.getCoefficient("A"), \
2718 self.getCoefficient("B"), \
2719 self.getCoefficient("C"), \
2720 self.getCoefficient("D"), \
2721 escript.Data(), \
2722 escript.Data(), \
2723 self.getCoefficient("d"), \
2724 escript.Data(),\
2725 self.getCoefficient("d_contact"), \
2726 escript.Data())
2727 self.getDomain().addPDEToSystem(operator,escript.Data(), \
2728 self.getCoefficient("A_reduced"), \
2729 self.getCoefficient("B_reduced"), \
2730 self.getCoefficient("C_reduced"), \
2731 self.getCoefficient("D_reduced"), \
2732 escript.Data(), \
2733 escript.Data(), \
2734 self.getCoefficient("d_reduced"), \
2735 escript.Data(),\
2736 self.getCoefficient("d_contact_reduced"), \
2737 escript.Data())
2738 self.insertConstraint(rhs_only=False)
2739 self.trace("New operator has been built.")
2740 self.validOperator()
2741 return (self.getCurrentOperator(), self.getCurrentRightHandSide())
2742
2743 def insertConstraint(self, rhs_only=False):
2744 """
2745 Applies the constraints defined by q and r to the PDE.
2746
2747 @param rhs_only: if True only the right hand side is altered by the
2748 constraint
2749 @type rhs_only: C{bool}
2750 """
2751 q=self.getCoefficient("q")
2752 r=self.getCoefficient("r")
2753 righthandside=self.getCurrentRightHandSide()
2754 operator=self.getCurrentOperator()
2755
2756 if not q.isEmpty():
2757 if r.isEmpty():
2758 r_s=self.createSolution()
2759 else:
2760 r_s=r
2761 if not rhs_only and not operator.isEmpty():
2762 if self.isUsingLumping():
2763 operator.copyWithMask(escript.Data(1.,q.getShape(),q.getFunctionSpace()),q)
2764 else:
2765 row_q=escript.Data(q,self.getFunctionSpaceForEquation())
2766 col_q=escript.Data(q,self.getFunctionSpaceForSolution())
2767 u=self.createSolution()
2768 u.copyWithMask(r_s,col_q)
2769 righthandside-=operator*u
2770 operator.nullifyRowsAndCols(row_q,col_q,1.)
2771 righthandside.copyWithMask(r_s,q)
2772
2773 def setValue(self,**coefficients):
2774 """
2775 Sets new values to coefficients.
2776
2777 @param coefficients: new values assigned to coefficients
2778 @keyword A: value for coefficient C{A}
2779 @type A: any type that can be cast to a L{Data<escript.Data>} object on
2780 L{Function<escript.Function>}
2781 @keyword A_reduced: value for coefficient C{A_reduced}
2782 @type A_reduced: any type that can be cast to a L{Data<escript.Data>}
2783 object on L{ReducedFunction<escript.ReducedFunction>}
2784 @keyword B: value for coefficient C{B}
2785 @type B: any type that can be cast to a L{Data<escript.Data>} object on
2786 L{Function<escript.Function>}
2787 @keyword B_reduced: value for coefficient C{B_reduced}
2788 @type B_reduced: any type that can be cast to a L{Data<escript.Data>}
2789 object on L{ReducedFunction<escript.ReducedFunction>}
2790 @keyword C: value for coefficient C{C}
2791 @type C: any type that can be cast to a L{Data<escript.Data>} object on
2792 L{Function<escript.Function>}
2793 @keyword C_reduced: value for coefficient C{C_reduced}
2794 @type C_reduced: any type that can be cast to a L{Data<escript.Data>}
2795 object on L{ReducedFunction<escript.ReducedFunction>}
2796 @keyword D: value for coefficient C{D}
2797 @type D: any type that can be cast to a L{Data<escript.Data>} object on
2798 L{Function<escript.Function>}
2799 @keyword D_reduced: value for coefficient C{D_reduced}
2800 @type D_reduced: any type that can be cast to a L{Data<escript.Data>}
2801 object on L{ReducedFunction<escript.ReducedFunction>}
2802 @keyword X: value for coefficient C{X}
2803 @type X: any type that can be cast to a L{Data<escript.Data>} object on
2804 L{Function<escript.Function>}
2805 @keyword X_reduced: value for coefficient C{X_reduced}
2806 @type X_reduced: any type that can be cast to a L{Data<escript.Data>}
2807 object on L{ReducedFunction<escript.ReducedFunction>}
2808 @keyword Y: value for coefficient C{Y}
2809 @type Y: any type that can be cast to a L{Data<escript.Data>} object on
2810 L{Function<escript.Function>}
2811 @keyword Y_reduced: value for coefficient C{Y_reduced}
2812 @type Y_reduced: any type that can be cast to a L{Data<escript.Data>}
2813 object on L{ReducedFunction<escript.Function>}
2814 @keyword d: value for coefficient C{d}
2815 @type d: any type that can be cast to a L{Data<escript.Data>} object on
2816 L{FunctionOnBoundary<escript.FunctionOnBoundary>}
2817 @keyword d_reduced: value for coefficient C{d_reduced}
2818 @type d_reduced: any type that can be cast to a L{Data<escript.Data>}
2819 object on L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}
2820 @keyword y: value for coefficient C{y}
2821 @type y: any type that can be cast to a L{Data<escript.Data>} object on
2822 L{FunctionOnBoundary<escript.FunctionOnBoundary>}
2823 @keyword d_contact: value for coefficient C{d_contact}
2824 @type d_contact: any type that can be cast to a L{Data<escript.Data>}
2825 object on L{FunctionOnContactOne<escript.FunctionOnContactOne>}
2826 or L{FunctionOnContactZero<escript.FunctionOnContactZero>}
2827 @keyword d_contact_reduced: value for coefficient C{d_contact_reduced}
2828 @type d_contact_reduced: any type that can be cast to a L{Data<escript.Data>}
2829 object on L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}
2830 or L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>}
2831 @keyword y_contact: value for coefficient C{y_contact}
2832 @type y_contact: any type that can be cast to a L{Data<escript.Data>}
2833 object on L{FunctionOnContactOne<escript.FunctionOnContactOne>}
2834 or L{FunctionOnContactZero<escript.FunctionOnContactZero>}
2835 @keyword y_contact_reduced: value for coefficient C{y_contact_reduced}
2836 @type y_contact_reduced: any type that can be cast to a L{Data<escript.Data>}
2837 object on L{ReducedFunctionOnContactOne<escript.FunctionOnContactOne>}
2838 or L{ReducedFunctionOnContactZero<escript.FunctionOnContactZero>}
2839 @keyword r: values prescribed to the solution at the locations of
2840 constraints
2841 @type r: any type that can be cast to a L{Data<escript.Data>} object on
2842 L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2843 depending on whether reduced order is used for the solution
2844 @keyword q: mask for location of constraints
2845 @type q: any type that can be cast to a L{Data<escript.Data>} object on
2846 L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2847 depending on whether reduced order is used for the
2848 representation of the equation
2849 @raise IllegalCoefficient: if an unknown coefficient keyword is used
2850 """
2851 super(LinearPDE,self).setValue(**coefficients)
2852 # check if the systrem is inhomogeneous:
2853 if len(coefficients)>0 and not self.isUsingLumping():
2854 q=self.getCoefficient("q")
2855 r=self.getCoefficient("r")
2856 if not q.isEmpty() and not r.isEmpty():
2857 if util.Lsup(q*r)>0.:
2858 self.trace("Inhomogeneous constraint detected.")
2859 self.invalidateSystem()
2860
2861
2862 def getResidual(self,u=None):
2863 """
2864 Returns the residual of u or the current solution if u is not present.
2865
2866 @param u: argument in the residual calculation. It must be representable
2867 in L{self.getFunctionSpaceForSolution()}. If u is not present
2868 or equals C{None} the current solution is used.
2869 @type u: L{Data<escript.Data>} or None
2870 @return: residual of u
2871 @rtype: L{Data<escript.Data>}
2872 """
2873 if u==None:
2874 return self.getOperator()*self.getSolution()-self.getRightHandSide()
2875 else:
2876 return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())-self.getRightHandSide()
2877
2878 def getFlux(self,u=None):
2879 """
2880 Returns the flux M{J} for a given M{u}.
2881
2882 M{J[i,j]=(A[i,j,k,l]+A_reduced[A[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]}
2883
2884 or
2885
2886 M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[l]+(B[j]+B_reduced[j])u-X[j]-X_reduced[j]}
2887
2888 @param u: argument in the flux. If u is not present or equals L{None} the
2889 current solution is used.
2890 @type u: L{Data<escript.Data>} or None
2891 @return: flux
2892 @rtype: L{Data<escript.Data>}
2893 """
2894 if u==None: u=self.getSolution()
2895 return util.tensormult(self.getCoefficient("A"),util.grad(u,Funtion(self.getDomain))) \
2896 +util.matrixmult(self.getCoefficient("B"),u) \
2897 -util.self.getCoefficient("X") \
2898 +util.tensormult(self.getCoefficient("A_reduced"),util.grad(u,ReducedFuntion(self.getDomain))) \
2899 +util.matrixmult(self.getCoefficient("B_reduced"),u) \
2900 -util.self.getCoefficient("X_reduced")
2901
2902
2903 class Poisson(LinearPDE):
2904 """
2905 Class to define a Poisson equation problem. This is generally a
2906 L{LinearPDE} of the form
2907
2908 M{-grad(grad(u)[j])[j] = f}
2909
2910 with natural boundary conditions
2911
2912 M{n[j]*grad(u)[j] = 0 }
2913
2914 and constraints:
2915
2916 M{u=0} where M{q>0}
2917
2918 """
2919
2920 def __init__(self,domain,debug=False):
2921 """
2922 Initializes a new Poisson equation.
2923
2924 @param domain: domain of the PDE
2925 @type domain: L{Domain<escript.Domain>}
2926 @param debug: if True debug information is printed
2927
2928 """
2929 super(Poisson, self).__init__(domain,1,1,debug)
2930 self.introduceCoefficients(
2931 f=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2932 f_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE))
2933 self.setSymmetryOn()
2934
2935 def setValue(self,**coefficients):
2936 """
2937 Sets new values to coefficients.
2938
2939 @param coefficients: new values assigned to coefficients
2940 @keyword f: value for right hand side M{f}
2941 @type f: any type that can be cast to a L{Scalar<escript.Scalar>} object
2942 on L{Function<escript.Function>}
2943 @keyword q: mask for location of constraints
2944 @type q: any type that can be cast to a rank zero L{Data<escript.Data>}
2945 object on L{Solution<escript.Solution>} or
2946 L{ReducedSolution<escript.ReducedSolution>} depending on whether
2947 reduced order is used for the representation of the equation
2948 @raise IllegalCoefficient: if an unknown coefficient keyword is used
2949 """
2950 super(Poisson, self).setValue(**coefficients)
2951
2952
2953 def getCoefficient(self,name):
2954 """
2955 Returns the value of the coefficient C{name} of the general PDE.
2956
2957 @param name: name of the coefficient requested
2958 @type name: C{string}
2959 @return: the value of the coefficient C{name}
2960 @rtype: L{Data<escript.Data>}
2961 @raise IllegalCoefficient: invalid coefficient name
2962 @note: This method is called by the assembling routine to map the Poisson
2963 equation onto the general PDE.
2964 """
2965 if name == "A" :
2966 return escript.Data(util.kronecker(self.getDim()),escript.Function(self.getDomain()))
2967 elif name == "Y" :
2968 return self.getCoefficient("f")
2969 elif name == "Y_reduced" :
2970 return self.getCoefficient("f_reduced")
2971 else:
2972 return super(Poisson, self).getCoefficient(name)
2973
2974 class Helmholtz(LinearPDE):
2975 """
2976 Class to define a Helmholtz equation problem. This is generally a
2977 L{LinearPDE} of the form
2978
2979 M{S{omega}*u - grad(k*grad(u)[j])[j] = f}
2980
2981 with natural boundary conditions
2982
2983 M{k*n[j]*grad(u)[j] = g- S{alpha}u }
2984
2985 and constraints:
2986
2987 M{u=r} where M{q>0}
2988
2989 """
2990
2991 def __init__(self,domain,debug=False):
2992 """
2993 Initializes a new Helmholtz equation.
2994
2995 @param domain: domain of the PDE
2996 @type domain: L{Domain<escript.Domain>}
2997 @param debug: if True debug information is printed
2998
2999 """
3000 super(Helmholtz, self).__init__(domain,1,1,debug)
3001 self.introduceCoefficients(
3002 omega=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.OPERATOR),
3003 k=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.OPERATOR),
3004 f=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
3005 f_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
3006 alpha=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.OPERATOR),
3007 g=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
3008 g_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE))
3009 self.setSymmetryOn()
3010
3011 def setValue(self,**coefficients):
3012 """
3013 Sets new values to coefficients.
3014
3015 @param coefficients: new values assigned to coefficients
3016 @keyword omega: value for coefficient M{S{omega}}
3017 @type omega: any type that can be cast to a L{Scalar<escript.Scalar>}
3018 object on L{Function<escript.Function>}
3019 @keyword k: value for coefficient M{k}
3020 @type k: any type that can be cast to a L{Scalar<escript.Scalar>} object
3021 on L{Function<escript.Function>}
3022 @keyword f: value for right hand side M{f}
3023 @type f: any type that can be cast to a L{Scalar<escript.Scalar>} object
3024 on L{Function<escript.Function>}
3025 @keyword alpha: value for right hand side M{S{alpha}}
3026 @type alpha: any type that can be cast to a L{Scalar<escript.Scalar>}
3027 object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}
3028 @keyword g: value for right hand side M{g}
3029 @type g: any type that can be cast to a L{Scalar<escript.Scalar>} object
3030 on L{FunctionOnBoundary<escript.FunctionOnBoundary>}
3031 @keyword r: prescribed values M{r} for the solution in constraints
3032 @type r: any type that can be cast to a L{Scalar<escript.Scalar>} object
3033 on L{Solution<escript.Solution>} or
3034 L{ReducedSolution<escript.ReducedSolution>} depending on whether
3035 reduced order is used for the representation of the equation
3036 @keyword q: mask for the location of constraints
3037 @type q: any type that can be cast to a L{Scalar<escript.Scalar>} object
3038 on L{Solution<escript.Solution>} or
3039 L{ReducedSolution<escript.ReducedSolution>} depending on whether
3040 reduced order is used for the representation of the equation
3041 @raise IllegalCoefficient: if an unknown coefficient keyword is used
3042 """
3043 super(Helmholtz, self).setValue(**coefficients)
3044
3045 def getCoefficient(self,name):
3046 """
3047 Returns the value of the coefficient C{name} of the general PDE.
3048
3049 @param name: name of the coefficient requested
3050 @type name: C{string}
3051 @return: the value of the coefficient C{name}
3052 @rtype: L{Data<escript.Data>}
3053 @raise IllegalCoefficient: invalid name
3054 """
3055 if name == "A" :
3056 if self.getCoefficient("k").isEmpty():
3057 return escript.Data(numpy.identity(self.getDim()),escript.Function(self.getDomain()))
3058 else:
3059 return escript.Data(numpy.identity(self.getDim()),escript.Function(self.getDomain()))*self.getCoefficient("k")
3060 elif name == "D" :
3061 return self.getCoefficient("omega")
3062 elif name == "Y" :
3063 return self.getCoefficient("f")
3064 elif name == "d" :
3065 return self.getCoefficient("alpha")
3066 elif name == "y" :
3067 return self.getCoefficient("g")
3068 elif name == "Y_reduced" :
3069 return self.getCoefficient("f_reduced")
3070 elif name == "y_reduced" :
3071 return self.getCoefficient("g_reduced")
3072 else:
3073 return super(Helmholtz, self).getCoefficient(name)
3074
3075 class LameEquation(LinearPDE):
3076 """
3077 Class to define a Lame equation problem. This problem is defined as:
3078
3079 M{-grad(S{mu}*(grad(u[i])[j]+grad(u[j])[i]))[j] - grad(S{lambda}*grad(u[k])[k])[j] = F_i -grad(S{sigma}[ij])[j] }
3080
3081 with natural boundary conditions:
3082
3083 M{n[j]*(S{mu}*(grad(u[i])[j]+grad(u[j])[i]) + S{lambda}*grad(u[k])[k]) = f_i +n[j]*S{sigma}[ij] }
3084
3085 and constraints:
3086
3087 M{u[i]=r[i]} where M{q[i]>0}
3088
3089 """
3090
3091 def __init__(self,domain,debug=False):
3092 """
3093 Initializes a new Lame equation.
3094
3095 @param domain: domain of the PDE
3096 @type domain: L{Domain<escript.Domain>}
3097 @param debug: if True debug information is printed
3098
3099 """
3100 super(LameEquation, self).__init__(domain,\
3101 domain.getDim(),domain.getDim(),debug)
3102 self.introduceCoefficients(lame_lambda=PDECoef(PDECoef.INTERIOR,(),PDECoef.OPERATOR),
3103 lame_mu=PDECoef(PDECoef.INTERIOR,(),PDECoef.OPERATOR),
3104 F=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
3105 sigma=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
3106 f=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE))
3107 self.setSymmetryOn()
3108
3109 def setValues(self,**coefficients):
3110 """
3111 Sets new values to coefficients.
3112
3113 @param coefficients: new values assigned to coefficients
3114 @keyword lame_mu: value for coefficient M{S{mu}}
3115 @type lame_mu: any type that can be cast to a L{Scalar<escript.Scalar>}
3116 object on L{Function<escript.Function>}
3117 @keyword lame_lambda: value for coefficient M{S{lambda}}
3118 @type lame_lambda: any type that can be cast to a L{Scalar<escript.Scalar>}
3119 object on L{Function<escript.Function>}
3120 @keyword F: value for internal force M{F}
3121 @type F: any type that can be cast to a L{Vector<escript.Vector>} object
3122 on L{Function<escript.Function>}
3123 @keyword sigma: value for initial stress M{S{sigma}}
3124 @type sigma: any type that can be cast to a L{Tensor<escript.Tensor>}
3125 object on L{Function<escript.Function>}
3126 @keyword f: value for external force M{f}
3127 @type f: any type that can be cast to a L{Vector<escript.Vector>} object
3128 on L{FunctionOnBoundary<escript.FunctionOnBoundary>}
3129 @keyword r: prescribed values M{r} for the solution in constraints
3130 @type r: any type that can be cast to a L{Vector<escript.Vector>} object
3131 on L{Solution<escript.Solution>} or
3132 L{ReducedSolution<escript.ReducedSolution>} depending on whether
3133 reduced order is used for the representation of the equation
3134 @keyword q: mask for the location of constraints
3135 @type q: any type that can be cast to a L{Vector<escript.Vector>} object
3136 on L{Solution<escript.Solution>} or
3137 L{ReducedSolution<escript.ReducedSolution>} depending on whether
3138 reduced order is used for the representation of the equation
3139 @raise IllegalCoefficient: if an unknown coefficient keyword is used
3140 """
3141 super(LameEquation, self).setValues(**coefficients)
3142
3143 def getCoefficient(self,name):
3144 """
3145 Returns the value of the coefficient C{name} of the general PDE.
3146
3147 @param name: name of the coefficient requested
3148 @type name: C{string}
3149 @return: the value of the coefficient C{name}
3150 @rtype: L{Data<escript.Data>}
3151 @raise IllegalCoefficient: invalid coefficient name
3152 """
3153 out =self.createCoefficient("A")
3154 if name == "A" :
3155 if self.getCoefficient("lame_lambda").isEmpty():
3156 if self.getCoefficient("lame_mu").isEmpty():
3157 pass
3158 else:
3159 for i in range(self.getDim()):
3160 for j in range(self.getDim()):
3161 out[i,j,j,i] += self.getCoefficient("lame_mu")
3162 out[i,j,i,j] += self.getCoefficient("lame_mu")
3163 else:
3164 if self.getCoefficient("lame_mu").isEmpty():
3165 for i in range(self.getDim()):
3166 for j in range(self.getDim()):
3167 out[i,i,j,j] += self.getCoefficient("lame_lambda")
3168 else:
3169 for i in range(self.getDim()):
3170 for j in range(self.getDim()):
3171 out[i,i,j,j] += self.getCoefficient("lame_lambda")
3172 out[i,j,j,i] += self.getCoefficient("lame_mu")
3173 out[i,j,i,j] += self.getCoefficient("lame_mu")
3174 return out
3175 elif name == "X" :
3176 return self.getCoefficient("sigma")
3177 elif name == "Y" :
3178 return self.getCoefficient("F")
3179 elif name == "y" :
3180 return self.getCoefficient("f")
3181 else:
3182 return super(LameEquation, self).getCoefficient(name)
3183
3184 def LinearSinglePDE(domain,debug=False):
3185 """
3186 Defines a single linear PDE.
3187
3188 @param domain: domain of the PDE
3189 @type domain: L{Domain<escript.Domain>}
3190 @param debug: if True debug information is printed
3191 @rtype: L{LinearPDE}
3192 """
3193 return LinearPDE(domain,numEquations=1,numSolutions=1,debug=debug)
3194
3195 def LinearPDESystem(domain,debug=False):
3196 """
3197 Defines a system of linear PDEs.
3198
3199 @param domain: domain of the PDEs
3200 @type domain: L{Domain<escript.Domain>}
3201 @param debug: if True debug information is printed
3202 @rtype: L{LinearPDE}
3203 """
3204 return LinearPDE(domain,numEquations=domain.getDim(),numSolutions=domain.getDim(),debug=debug)
3205
3206
3207 class TransportPDE(LinearProblem):
3208 """
3209 This class is used to define a transport problem given by a general linear,
3210 time dependent, second order PDE for an unknown, non-negative,
3211 time-dependent function M{u} on a given domain defined through a
3212 L{Domain<escript.Domain>} object.
3213
3214 For a single equation with a solution with a single component the transport
3215 problem is defined in the following form:
3216
3217 M{(M+M_reduced)*u_t=-(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)}
3218
3219 where M{u_t} denotes the time derivative of M{u} and M{grad(F)} denotes the
3220 spatial derivative of M{F}. Einstein's summation convention, ie. summation
3221 over indexes appearing twice in a term of a sum performed, is used.
3222 The coefficients M{M}, M{A}, M{B}, M{C}, M{D}, M{X} and M{Y} have to be
3223 specified through L{Data<escript.Data>} objects in L{Function<escript.Function>}
3224 and the coefficients M{M_reduced}, M{A_reduced}, M{B_reduced}, M{C_reduced},
3225 M{D_reduced}, M{X_reduced} and M{Y_reduced} have to be specified through
3226 L{Data<escript.Data>} objects in L{ReducedFunction<escript.ReducedFunction>}.
3227 It is also allowed to use objects that can be converted into such
3228 L{Data<escript.Data>} objects. M{M} and M{M_reduced} are scalar, M{A} and
3229 M{A_reduced} are rank two, M{B}, M{C}, M{X}, M{B_reduced}, M{C_reduced} and
3230 M{X_reduced} are rank one and M{D}, M{D_reduced}, M{Y} and M{Y_reduced} are
3231 scalar.
3232
3233 The following natural boundary conditions are considered:
3234
3235 M{n[j]*((A[i,j]+A_reduced[i,j])*grad(u)[l]+(B+B_reduced)[j]*u+X[j]+X_reduced[j])+(d+d_reduced)*u+y+y_reduced=(m+m_reduced)*u_t}
3236
3237 where M{n} is the outer normal field. Notice that the coefficients M{A},
3238 M{A_reduced}, M{B}, M{B_reduced}, M{X} and M{X_reduced} are defined in the
3239 transport problem. The coefficients M{m}, M{d} and M{y} are each a scalar in
3240 L{FunctionOnBoundary<escript.FunctionOnBoundary>} and the coefficients
3241 M{m_reduced}, M{d_reduced} and M{y_reduced} are each a scalar in
3242 L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
3243
3244 Constraints for the solution prescribing the value of the solution at
3245 certain locations in the domain have the form
3246
3247 M{u_t=r} where M{q>0}
3248
3249 M{r} and M{q} are each scalar where M{q} is the characteristic function
3250 defining where the constraint is applied. The constraints override any other
3251 condition set by the transport problem or the boundary condition.
3252
3253 The transport problem is symmetrical if
3254
3255 M{A[i,j]=A[j,i]} and M{B[j]=C[j]} and M{A_reduced[i,j]=A_reduced[j,i]} and
3256 M{B_reduced[j]=C_reduced[j]}
3257
3258 For a system and a solution with several components the transport problem
3259 has the form
3260
3261 M{(M[i,k]+M_reduced[i,k])*u[k]_t=-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] }
3262
3263 M{A} and M{A_reduced} are of rank four, M{B}, M{B_reduced}, M{C} and
3264 M{C_reduced} are each of rank three, M{M}, M{M_reduced}, M{D}, M{D_reduced},
3265 M{X_reduced} and M{X} are each of rank two and M{Y} and M{Y_reduced} are of
3266 rank one. The natural boundary conditions take the form:
3267
3268 M{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]+X[i,j]+X_reduced[i,j])+(d[i,k]+d_reduced[i,k])*u[k]+y[i]+y_reduced[i]= (m[i,k]+m_reduced[i,k])*u[k]_t}
3269
3270 The coefficient M{d} and M{m} are of rank two and M{y} is of rank one with
3271 all in L{FunctionOnBoundary<escript.FunctionOnBoundary>}. The coefficients
3272 M{d_reduced} and M{m_reduced} are of rank two and M{y_reduced} is of rank
3273 one all in L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
3274
3275 Constraints take the form
3276
3277 M{u[i]_t=r[i]} where M{q[i]>0}
3278
3279 M{r} and M{q} are each rank one. Notice that at some locations not
3280 necessarily all components must have a constraint.
3281
3282 The transport problem is symmetrical if
3283
3284 - M{M[i,k]=M[i,k]}
3285 - M{M_reduced[i,k]=M_reduced[i,k]}
3286 - M{A[i,j,k,l]=A[k,l,i,j]}
3287 - M{A_reduced[i,j,k,l]=A_reduced[k,l,i,j]}
3288 - M{B[i,j,k]=C[k,i,j]}
3289 - M{B_reduced[i,j,k]=C_reduced[k,i,j]}
3290 - M{D[i,k]=D[i,k]}
3291 - M{D_reduced[i,k]=D_reduced[i,k]}
3292 - M{m[i,k]=m[k,i]}
3293 - M{m_reduced[i,k]=m_reduced[k,i]}
3294 - M{d[i,k]=d[k,i]}
3295 - M{d_reduced[i,k]=d_reduced[k,i]}
3296
3297 L{TransportPDE} also supports solution discontinuities over a contact region
3298 in the domain. To specify the conditions across the discontinuity we are
3299 using the generalised flux M{J} which, in the case of a system of PDEs and
3300 several components of the solution, is defined as
3301
3302 M{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]}
3303
3304 For the case of single solution component and single PDE M{J} is defined as
3305
3306 M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[j]+(B[i]+B_reduced[i])*u+X[i]+X_reduced[i]}
3307
3308 In the context of discontinuities M{n} denotes the normal on the
3309 discontinuity pointing from side 0 towards side 1 calculated from
3310 L{getNormal<escript.FunctionSpace.getNormal>} of L{FunctionOnContactZero<escript.FunctionOnContactZero>}.
3311 For a system of transport problems the contact condition takes the form
3312
3313 M{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]}
3314
3315 where M{J0} and M{J1} are the fluxes on side 0 and side 1 of the
3316 discontinuity, respectively. M{jump(u)}, which is the difference of the
3317 solution at side 1 and at side 0, denotes the jump of M{u} across
3318 discontinuity along the normal calculated by L{jump<util.jump>}.
3319 The coefficient M{d_contact} is of rank two and M{y_contact} is of rank one
3320 both in L{FunctionOnContactZero<escript.FunctionOnContactZero>} or L{FunctionOnContactOne<escript.FunctionOnContactOne>}.
3321 The coefficient M{d_contact_reduced} is of rank two and M{y_contact_reduced}
3322 is of rank one both in L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}.
3323 In case of a single PDE and a single component solution the contact
3324 condition takes the form
3325
3326 M{n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)}
3327
3328 In this case the coefficient M{d_contact} and M{y_contact} are each scalar
3329 both in L{FunctionOnContactZero<escript.FunctionOnContactZero>} or
3330 L{FunctionOnContactOne<escript.FunctionOnContactOne>} and the coefficient
3331 M{d_contact_reduced} and M{y_contact_reduced} are each scalar both in
3332 L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or
3333 L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}.
3334
3335 Typical usage::
3336
3337 p = TransportPDE(dom)
3338 p.setValue(M=1., C=[-1.,0.])
3339 p.setInitialSolution(u=exp(-length(dom.getX()-[0.1,0.1])**2)
3340 t = 0
3341 dt = 0.1
3342 while (t < 1.):
3343 u = p.solve(dt)
3344
3345 """
3346 def __init__(self,domain,numEquations=None,numSolutions=None, useBackwardEuler=False, debug=False):
3347 """
3348 Initializes a transport problem.
3349
3350 @param domain: domain of the PDE
3351 @type domain: L{Domain<escript.Domain>}
3352 @param numEquations: number of equations. If C{None} the number of
3353 equations is extracted from the coefficients.
3354 @param numSolutions: number of solution components. If C{None} the number
3355 of solution components is extracted from the
3356 coefficients.
3357 @param debug: if True debug information is printed
3358 @param useBackwardEuler: if set the backward Euler scheme is used. Otherwise the Crank-Nicholson scheme is applied. Note that backward Euler scheme will return a safe time step size which is practically infinity as the scheme is unconditional unstable. The Crank-Nicholson scheme provides a higher accuracy but requires to limit the time step size to be stable.
3359 @type useBackwardEuler: C{bool}
3360 """
3361 if useBackwardEuler:
3362 self.__useBackwardEuler=True
3363 else:
3364 self.__useBackwardEuler=False
3365 super(TransportPDE, self).__init__(domain,numEquations,numSolutions,debug)
3366
3367 self.setConstraintWeightingFactor()
3368 #
3369 # the coefficients of the transport problem
3370 #
3371 self.introduceCoefficients(
3372 M=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
3373 A=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
3374 B=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
3375 C=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
3376 D=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
3377 X=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
3378 Y=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
3379 m=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
3380 d=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
3381 y=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
3382 d_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
3383 y_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
3384 M_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
3385 A_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
3386 B_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
3387 C_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
3388 D_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
3389 X_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
3390 Y_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
3391 m_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
3392 d_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
3393 y_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
3394 d_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
3395 y_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
3396 r=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.RIGHTHANDSIDE),
3397 q=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.BOTH) )
3398
3399 def __str__(self):
3400 """
3401 Returns the string representation of the transport problem.
3402
3403 @return: a simple representation of the transport problem
3404 @rtype: C{str}
3405 """
3406 return "<TransportPDE %d>"%id(self)
3407
3408 def useBackwardEu