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