/[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 2551 - (show annotations)
Thu Jul 23 09:19:15 2009 UTC (10 years, 1 month ago) by gross
File MIME type: text/x-python
File size: 153768 byte(s)
a problem with the sparse matrix unrolling fixed.
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 if self.getOperatorType() == self.getRequiredOperatorType():
2068 self.__operator.resetValues()
2069 else:
2070 self.__operator=self.createOperator()
2071 self.__operator_type=self.getRequiredOperatorType()
2072 self.trace("Operator reset to zero")
2073
2074 def getCurrentOperator(self):
2075 """
2076 Returns the operator in its current state.
2077 """
2078 return self.__operator
2079
2080 def setValue(self,**coefficients):
2081 """
2082 Sets new values to coefficients.
2083
2084 @raise IllegalCoefficient: if an unknown coefficient keyword is used
2085 """
2086 # check if the coefficients are legal:
2087 for i in coefficients.iterkeys():
2088 if not self.hasCoefficient(i):
2089 raise IllegalCoefficient,"Attempt to set unknown coefficient %s"%i
2090 # if the number of unknowns or equations is still unknown we try to estimate them:
2091 if self.__numEquations==None or self.__numSolutions==None:
2092 for i,d in coefficients.iteritems():
2093 if hasattr(d,"shape"):
2094 s=d.shape
2095 elif hasattr(d,"getShape"):
2096 s=d.getShape()
2097 else:
2098 s=numpy.array(d).shape
2099 if s!=None:
2100 # get number of equations and number of unknowns:
2101 res=self.__COEFFICIENTS[i].estimateNumEquationsAndNumSolutions(self.getDomain(),s)
2102 if res==None:
2103 raise IllegalCoefficientValue,"Illegal shape %s of coefficient %s"%(s,i)
2104 else:
2105 if self.__numEquations==None: self.__numEquations=res[0]
2106 if self.__numSolutions==None: self.__numSolutions=res[1]
2107 if self.__numEquations==None: raise UndefinedPDEError,"unidentified number of equations"
2108 if self.__numSolutions==None: raise UndefinedPDEError,"unidentified number of solutions"
2109 # now we check the shape of the coefficient if numEquations and numSolutions are set:
2110 for i,d in coefficients.iteritems():
2111 try:
2112 self.__COEFFICIENTS[i].setValue(self.getDomain(),
2113 self.getNumEquations(),self.getNumSolutions(),
2114 self.reduceEquationOrder(),self.reduceSolutionOrder(),d)
2115 self.alteredCoefficient(i)
2116 except IllegalCoefficientFunctionSpace,m:
2117 # if the function space is wrong then we try the reduced version:
2118 i_red=i+"_reduced"
2119 if (not i_red in coefficients.keys()) and i_red in self.__COEFFICIENTS.keys():
2120 try:
2121 self.__COEFFICIENTS[i_red].setValue(self.getDomain(),
2122 self.getNumEquations(),self.getNumSolutions(),
2123 self.reduceEquationOrder(),self.reduceSolutionOrder(),d)
2124 self.alteredCoefficient(i_red)
2125 except IllegalCoefficientValue,m:
2126 raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
2127 except IllegalCoefficientFunctionSpace,m:
2128 raise IllegalCoefficientFunctionSpace("Coefficient %s:%s"%(i,m))
2129 else:
2130 raise IllegalCoefficientFunctionSpace("Coefficient %s:%s"%(i,m))
2131 except IllegalCoefficientValue,m:
2132 raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
2133 self.__altered_coefficients=True
2134
2135 # ==========================================================================
2136 # methods that are typically overwritten when implementing a particular
2137 # linear problem
2138 # ==========================================================================
2139 def getRequiredOperatorType(self):
2140 """
2141 Returns the system type which needs to be used by the current set up.
2142
2143 @note: Typically this method is overwritten when implementing a
2144 particular linear problem.
2145 """
2146 return None
2147
2148 def createOperator(self):
2149 """
2150 Returns an instance of a new operator.
2151
2152 @note: This method is overwritten when implementing a particular
2153 linear problem.
2154 """
2155 return escript.Operator()
2156
2157 def checkSymmetry(self,verbose=True):
2158 """
2159 Tests the PDE for symmetry.
2160
2161 @param verbose: if set to True or not present a report on coefficients
2162 which break the symmetry is printed
2163 @type verbose: C{bool}
2164 @return: True if the problem is symmetric
2165 @rtype: C{bool}
2166 @note: Typically this method is overwritten when implementing a
2167 particular linear problem.
2168 """
2169 out=True
2170 return out
2171
2172 def getSolution(self,**options):
2173 """
2174 Returns the solution of the problem.
2175
2176 @return: the solution
2177 @rtype: L{Data<escript.Data>}
2178
2179 @note: This method is overwritten when implementing a particular
2180 linear problem.
2181 """
2182 return self.getCurrentSolution()
2183
2184 def getSystem(self):
2185 """
2186 Returns the operator and right hand side of the PDE.
2187
2188 @return: the discrete version of the PDE
2189 @rtype: C{tuple} of L{Operator,<escript.Operator>} and L{Data<escript.Data>}.
2190
2191 @note: This method is overwritten when implementing a particular
2192 linear problem.
2193 """
2194 return (self.getCurrentOperator(), self.getCurrentRightHandSide())
2195
2196 class LinearPDE(LinearProblem):
2197 """
2198 This class is used to define a general linear, steady, second order PDE
2199 for an unknown function M{u} on a given domain defined through a
2200 L{Domain<escript.Domain>} object.
2201
2202 For a single PDE having a solution with a single component the linear PDE
2203 is defined in the following form:
2204
2205 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)}
2206
2207 where M{grad(F)} denotes the spatial derivative of M{F}. Einstein's
2208 summation convention, ie. summation over indexes appearing twice in a term
2209 of a sum performed, is used.
2210 The coefficients M{A}, M{B}, M{C}, M{D}, M{X} and M{Y} have to be specified
2211 through L{Data<escript.Data>} objects in L{Function<escript.Function>} and
2212 the coefficients M{A_reduced}, M{B_reduced}, M{C_reduced}, M{D_reduced},
2213 M{X_reduced} and M{Y_reduced} have to be specified through
2214 L{Data<escript.Data>} objects in L{ReducedFunction<escript.ReducedFunction>}.
2215 It is also allowed to use objects that can be converted into such
2216 L{Data<escript.Data>} objects. M{A} and M{A_reduced} are rank two, M{B},
2217 M{C}, M{X}, M{B_reduced}, M{C_reduced} and M{X_reduced} are rank one and
2218 M{D}, M{D_reduced}, M{Y} and M{Y_reduced} are scalar.
2219
2220 The following natural boundary conditions are considered:
2221
2222 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}
2223
2224 where M{n} is the outer normal field. Notice that the coefficients M{A},
2225 M{A_reduced}, M{B}, M{B_reduced}, M{X} and M{X_reduced} are defined in the
2226 PDE. The coefficients M{d} and M{y} are each a scalar in
2227 L{FunctionOnBoundary<escript.FunctionOnBoundary>} and the coefficients
2228 M{d_reduced} and M{y_reduced} are each a scalar in
2229 L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
2230
2231 Constraints for the solution prescribe the value of the solution at certain
2232 locations in the domain. They have the form
2233
2234 M{u=r} where M{q>0}
2235
2236 M{r} and M{q} are each scalar where M{q} is the characteristic function
2237 defining where the constraint is applied. The constraints override any
2238 other condition set by the PDE or the boundary condition.
2239
2240 The PDE is symmetrical if
2241
2242 M{A[i,j]=A[j,i]} and M{B[j]=C[j]} and M{A_reduced[i,j]=A_reduced[j,i]}
2243 and M{B_reduced[j]=C_reduced[j]}
2244
2245 For a system of PDEs and a solution with several components the PDE has the
2246 form
2247
2248 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] }
2249
2250 M{A} and M{A_reduced} are of rank four, M{B}, M{B_reduced}, M{C} and
2251 M{C_reduced} are each of rank three, M{D}, M{D_reduced}, M{X_reduced} and
2252 M{X} are each of rank two and M{Y} and M{Y_reduced} are of rank one.
2253 The natural boundary conditions take the form:
2254
2255 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]}
2256
2257 The coefficient M{d} is of rank two and M{y} is of rank one both in
2258 L{FunctionOnBoundary<escript.FunctionOnBoundary>}. The coefficients
2259 M{d_reduced} is of rank two and M{y_reduced} is of rank one both in
2260 L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
2261
2262 Constraints take the form
2263
2264 M{u[i]=r[i]} where M{q[i]>0}
2265
2266 M{r} and M{q} are each rank one. Notice that at some locations not
2267 necessarily all components must have a constraint.
2268
2269 The system of PDEs is symmetrical if
2270
2271 - M{A[i,j,k,l]=A[k,l,i,j]}
2272 - M{A_reduced[i,j,k,l]=A_reduced[k,l,i,j]}
2273 - M{B[i,j,k]=C[k,i,j]}
2274 - M{B_reduced[i,j,k]=C_reduced[k,i,j]}
2275 - M{D[i,k]=D[i,k]}
2276 - M{D_reduced[i,k]=D_reduced[i,k]}
2277 - M{d[i,k]=d[k,i]}
2278 - M{d_reduced[i,k]=d_reduced[k,i]}
2279
2280 L{LinearPDE} also supports solution discontinuities over a contact region
2281 in the domain. To specify the conditions across the discontinuity we are
2282 using the generalised flux M{J} which, in the case of a system of PDEs
2283 and several components of the solution, is defined as
2284
2285 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]}
2286
2287 For the case of single solution component and single PDE M{J} is defined as
2288
2289 M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[j]+(B[i]+B_reduced[i])*u-X[i]-X_reduced[i]}
2290
2291 In the context of discontinuities M{n} denotes the normal on the
2292 discontinuity pointing from side 0 towards side 1 calculated from
2293 L{getNormal<escript.FunctionSpace.getNormal>} of L{FunctionOnContactZero<escript.FunctionOnContactZero>}.
2294 For a system of PDEs the contact condition takes the form
2295
2296 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]}
2297
2298 where M{J0} and M{J1} are the fluxes on side 0 and side 1 of the
2299 discontinuity, respectively. M{jump(u)}, which is the difference of the
2300 solution at side 1 and at side 0, denotes the jump of M{u} across
2301 discontinuity along the normal calculated by L{jump<util.jump>}.
2302 The coefficient M{d_contact} is of rank two and M{y_contact} is of rank one
2303 both in L{FunctionOnContactZero<escript.FunctionOnContactZero>} or
2304 L{FunctionOnContactOne<escript.FunctionOnContactOne>}.
2305 The coefficient M{d_contact_reduced} is of rank two and M{y_contact_reduced}
2306 is of rank one both in L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>}
2307 or L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}.
2308 In case of a single PDE and a single component solution the contact
2309 condition takes the form
2310
2311 M{n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)}
2312
2313 In this case the coefficient M{d_contact} and M{y_contact} are each scalar
2314 both in L{FunctionOnContactZero<escript.FunctionOnContactZero>} or
2315 L{FunctionOnContactOne<escript.FunctionOnContactOne>} and the coefficient
2316 M{d_contact_reduced} and M{y_contact_reduced} are each scalar both in
2317 L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or
2318 L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}.
2319
2320 Typical usage::
2321
2322 p = LinearPDE(dom)
2323 p.setValue(A=kronecker(dom), D=1, Y=0.5)
2324 u = p.getSolution()
2325
2326 """
2327
2328 def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):
2329 """
2330 Initializes a new linear PDE.
2331
2332 @param domain: domain of the PDE
2333 @type domain: L{Domain<escript.Domain>}
2334 @param numEquations: number of equations. If C{None} the number of
2335 equations is extracted from the PDE coefficients.
2336 @param numSolutions: number of solution components. If C{None} the number
2337 of solution components is extracted from the PDE
2338 coefficients.
2339 @param debug: if True debug information is printed
2340
2341 """
2342 super(LinearPDE, self).__init__(domain,numEquations,numSolutions,debug)
2343 #
2344 # the coefficients of the PDE:
2345 #
2346 self.introduceCoefficients(
2347 A=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2348 B=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2349 C=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2350 D=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2351 X=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
2352 Y=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2353 d=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2354 y=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2355 d_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2356 y_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2357 A_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2358 B_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2359 C_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2360 D_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2361 X_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
2362 Y_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2363 d_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2364 y_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2365 d_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2366 y_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2367 r=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.RIGHTHANDSIDE),
2368 q=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.BOTH) )
2369
2370 def __str__(self):
2371 """
2372 Returns the string representation of the PDE.
2373
2374 @return: a simple representation of the PDE
2375 @rtype: C{str}
2376 """
2377 return "<LinearPDE %d>"%id(self)
2378
2379 def getRequiredOperatorType(self):
2380 """
2381 Returns the system type which needs to be used by the current set up.
2382 """
2383 solver_options=self.getSolverOptions()
2384 return self.getDomain().getSystemMatrixTypeId(solver_options.getSolverMethod(), solver_options.getPreconditioner(),solver_options.getPackage(), solver_options.isSymmetric())
2385
2386 def checkSymmetry(self,verbose=True):
2387 """
2388 Tests the PDE for symmetry.
2389
2390 @param verbose: if set to True or not present a report on coefficients
2391 which break the symmetry is printed.
2392 @type verbose: C{bool}
2393 @return: True if the PDE is symmetric
2394 @rtype: L{bool}
2395 @note: This is a very expensive operation. It should be used for
2396 degugging only! The symmetry flag is not altered.
2397 """
2398 out=True
2399 out=out and self.checkSymmetricTensor("A", verbose)
2400 out=out and self.checkSymmetricTensor("A_reduced", verbose)
2401 out=out and self.checkReciprocalSymmetry("B","C", verbose)
2402 out=out and self.checkReciprocalSymmetry("B_reduced","C_reduced", verbose)
2403 out=out and self.checkSymmetricTensor("D", verbose)
2404 out=out and self.checkSymmetricTensor("D_reduced", verbose)
2405 out=out and self.checkSymmetricTensor("d", verbose)
2406 out=out and self.checkSymmetricTensor("d_reduced", verbose)
2407 out=out and self.checkSymmetricTensor("d_contact", verbose)
2408 out=out and self.checkSymmetricTensor("d_contact_reduced", verbose)
2409 return out
2410
2411 def createOperator(self):
2412 """
2413 Returns an instance of a new operator.
2414 """
2415 optype=self.getRequiredOperatorType()
2416 self.trace("New operator of type %s is allocated."%optype)
2417 return self.getDomain().newOperator( \
2418 self.getNumEquations(), \
2419 self.getFunctionSpaceForEquation(), \
2420 self.getNumSolutions(), \
2421 self.getFunctionSpaceForSolution(), \
2422 optype)
2423
2424 def getSolution(self):
2425 """
2426 Returns the solution of the PDE.
2427
2428 @return: the solution
2429 @rtype: L{Data<escript.Data>}
2430 """
2431 option_class=self.getSolverOptions()
2432 if not self.isSolutionValid():
2433 mat,f=self.getSystem()
2434 if self.isUsingLumping():
2435 self.setSolution(f*1/mat)
2436 else:
2437 self.trace("PDE is resolved.")
2438 self.trace("solver options: %s"%str(option_class))
2439 self.setSolution(mat.solve(f,option_class))
2440 return self.getCurrentSolution()
2441
2442 def getSystem(self):
2443 """
2444 Returns the operator and right hand side of the PDE.
2445
2446 @return: the discrete version of the PDE
2447 @rtype: C{tuple} of L{Operator,<escript.Operator>} and
2448 L{Data<escript.Data>}
2449 """
2450 if not self.isOperatorValid() or not self.isRightHandSideValid():
2451 if self.isUsingLumping():
2452 if not self.isOperatorValid():
2453 if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():
2454 raise TypeError,"Lumped matrix requires same order for equations and unknowns"
2455 if not self.getCoefficient("A").isEmpty():
2456 raise ValueError,"coefficient A in lumped matrix may not be present."
2457 if not self.getCoefficient("B").isEmpty():
2458 raise ValueError,"coefficient B in lumped matrix may not be present."
2459 if not self.getCoefficient("C").isEmpty():
2460 raise ValueError,"coefficient C in lumped matrix may not be present."
2461 if not self.getCoefficient("d_contact").isEmpty():
2462 raise ValueError,"coefficient d_contact in lumped matrix may not be present."
2463 if not self.getCoefficient("A_reduced").isEmpty():
2464 raise ValueError,"coefficient A_reduced in lumped matrix may not be present."
2465 if not self.getCoefficient("B_reduced").isEmpty():
2466 raise ValueError,"coefficient B_reduced in lumped matrix may not be present."
2467 if not self.getCoefficient("C_reduced").isEmpty():
2468 raise ValueError,"coefficient C_reduced in lumped matrix may not be present."
2469 if not self.getCoefficient("d_contact_reduced").isEmpty():
2470 raise ValueError,"coefficient d_contact_reduced in lumped matrix may not be present."
2471 D=self.getCoefficient("D")
2472 d=self.getCoefficient("d")
2473 D_reduced=self.getCoefficient("D_reduced")
2474 d_reduced=self.getCoefficient("d_reduced")
2475 if not D.isEmpty():
2476 if self.getNumSolutions()>1:
2477 D_times_e=util.matrix_mult(D,numpy.ones((self.getNumSolutions(),)))
2478 else:
2479 D_times_e=D
2480 else:
2481 D_times_e=escript.Data()
2482 if not d.isEmpty():
2483 if self.getNumSolutions()>1:
2484 d_times_e=util.matrix_mult(d,numpy.ones((self.getNumSolutions(),)))
2485 else:
2486 d_times_e=d
2487 else:
2488 d_times_e=escript.Data()
2489
2490 if not D_reduced.isEmpty():
2491 if self.getNumSolutions()>1:
2492 D_reduced_times_e=util.matrix_mult(D_reduced,numpy.ones((self.getNumSolutions(),)))
2493 else:
2494 D_reduced_times_e=D_reduced
2495 else:
2496 D_reduced_times_e=escript.Data()
2497 if not d_reduced.isEmpty():
2498 if self.getNumSolutions()>1:
2499 d_reduced_times_e=util.matrix_mult(d_reduced,numpy.ones((self.getNumSolutions(),)))
2500 else:
2501 d_reduced_times_e=d_reduced
2502 else:
2503 d_reduced_times_e=escript.Data()
2504
2505 self.resetOperator()
2506 operator=self.getCurrentOperator()
2507 if False and hasattr(self.getDomain(), "addPDEToLumpedSystem") :
2508 self.getDomain().addPDEToLumpedSystem(operator, D_times_e, d_times_e)
2509 self.getDomain().addPDEToLumpedSystem(operator, D_reduced_times_e, d_reduced_times_e)
2510 else:
2511 self.getDomain().addPDEToRHS(operator, \
2512 escript.Data(), \
2513 D_times_e, \
2514 d_times_e,\
2515 escript.Data())
2516 self.getDomain().addPDEToRHS(operator, \
2517 escript.Data(), \
2518 D_reduced_times_e, \
2519 d_reduced_times_e,\
2520 escript.Data())
2521 self.trace("New lumped operator has been built.")
2522 if not self.isRightHandSideValid():
2523 self.resetRightHandSide()
2524 righthandside=self.getCurrentRightHandSide()
2525 self.getDomain().addPDEToRHS(righthandside, \
2526 self.getCoefficient("X"), \
2527 self.getCoefficient("Y"),\
2528 self.getCoefficient("y"),\
2529 self.getCoefficient("y_contact"))
2530 self.getDomain().addPDEToRHS(righthandside, \
2531 self.getCoefficient("X_reduced"), \
2532 self.getCoefficient("Y_reduced"),\
2533 self.getCoefficient("y_reduced"),\
2534 self.getCoefficient("y_contact_reduced"))
2535 self.trace("New right hand side has been built.")
2536 self.validRightHandSide()
2537 self.insertConstraint(rhs_only=False)
2538 self.validOperator()
2539 else:
2540 if not self.isOperatorValid() and not self.isRightHandSideValid():
2541 self.resetRightHandSide()
2542 righthandside=self.getCurrentRightHandSide()
2543 self.resetOperator()
2544 operator=self.getCurrentOperator()
2545 self.getDomain().addPDEToSystem(operator,righthandside, \
2546 self.getCoefficient("A"), \
2547 self.getCoefficient("B"), \
2548 self.getCoefficient("C"), \
2549 self.getCoefficient("D"), \
2550 self.getCoefficient("X"), \
2551 self.getCoefficient("Y"), \
2552 self.getCoefficient("d"), \
2553 self.getCoefficient("y"), \
2554 self.getCoefficient("d_contact"), \
2555 self.getCoefficient("y_contact"))
2556 self.getDomain().addPDEToSystem(operator,righthandside, \
2557 self.getCoefficient("A_reduced"), \
2558 self.getCoefficient("B_reduced"), \
2559 self.getCoefficient("C_reduced"), \
2560 self.getCoefficient("D_reduced"), \
2561 self.getCoefficient("X_reduced"), \
2562 self.getCoefficient("Y_reduced"), \
2563 self.getCoefficient("d_reduced"), \
2564 self.getCoefficient("y_reduced"), \
2565 self.getCoefficient("d_contact_reduced"), \
2566 self.getCoefficient("y_contact_reduced"))
2567 self.insertConstraint(rhs_only=False)
2568 self.trace("New system has been built.")
2569 self.validOperator()
2570 self.validRightHandSide()
2571 elif not self.isRightHandSideValid():
2572 self.resetRightHandSide()
2573 righthandside=self.getCurrentRightHandSide()
2574 self.getDomain().addPDEToRHS(righthandside,
2575 self.getCoefficient("X"), \
2576 self.getCoefficient("Y"),\
2577 self.getCoefficient("y"),\
2578 self.getCoefficient("y_contact"))
2579 self.getDomain().addPDEToRHS(righthandside,
2580 self.getCoefficient("X_reduced"), \
2581 self.getCoefficient("Y_reduced"),\
2582 self.getCoefficient("y_reduced"),\
2583 self.getCoefficient("y_contact_reduced"))
2584 self.insertConstraint(rhs_only=True)
2585 self.trace("New right hand side has been built.")
2586 self.validRightHandSide()
2587 elif not self.isOperatorValid():
2588 self.resetOperator()
2589 operator=self.getCurrentOperator()
2590 self.getDomain().addPDEToSystem(operator,escript.Data(), \
2591 self.getCoefficient("A"), \
2592 self.getCoefficient("B"), \
2593 self.getCoefficient("C"), \
2594 self.getCoefficient("D"), \
2595 escript.Data(), \
2596 escript.Data(), \
2597 self.getCoefficient("d"), \
2598 escript.Data(),\
2599 self.getCoefficient("d_contact"), \
2600 escript.Data())
2601 self.getDomain().addPDEToSystem(operator,escript.Data(), \
2602 self.getCoefficient("A_reduced"), \
2603 self.getCoefficient("B_reduced"), \
2604 self.getCoefficient("C_reduced"), \
2605 self.getCoefficient("D_reduced"), \
2606 escript.Data(), \
2607 escript.Data(), \
2608 self.getCoefficient("d_reduced"), \
2609 escript.Data(),\
2610 self.getCoefficient("d_contact_reduced"), \
2611 escript.Data())
2612 self.insertConstraint(rhs_only=False)
2613 self.trace("New operator has been built.")
2614 self.validOperator()
2615 self.setSystemStatus()
2616 self.trace("System status is %s."%self.getSystemStatus())
2617 return (self.getCurrentOperator(), self.getCurrentRightHandSide())
2618
2619 def insertConstraint(self, rhs_only=False):
2620 """
2621 Applies the constraints defined by q and r to the PDE.
2622
2623 @param rhs_only: if True only the right hand side is altered by the
2624 constraint
2625 @type rhs_only: C{bool}
2626 """
2627 q=self.getCoefficient("q")
2628 r=self.getCoefficient("r")
2629 righthandside=self.getCurrentRightHandSide()
2630 operator=self.getCurrentOperator()
2631
2632 if not q.isEmpty():
2633 if r.isEmpty():
2634 r_s=self.createSolution()
2635 else:
2636 r_s=r
2637 if not rhs_only and not operator.isEmpty():
2638 if self.isUsingLumping():
2639 operator.copyWithMask(escript.Data(1.,q.getShape(),q.getFunctionSpace()),q)
2640 else:
2641 row_q=escript.Data(q,self.getFunctionSpaceForEquation())
2642 col_q=escript.Data(q,self.getFunctionSpaceForSolution())
2643 u=self.createSolution()
2644 u.copyWithMask(r_s,col_q)
2645 righthandside-=operator*u
2646 operator.nullifyRowsAndCols(row_q,col_q,1.)
2647 righthandside.copyWithMask(r_s,q)
2648
2649 def setValue(self,**coefficients):
2650 """
2651 Sets new values to coefficients.
2652
2653 @param coefficients: new values assigned to coefficients
2654 @keyword A: value for coefficient C{A}
2655 @type A: any type that can be cast to a L{Data<escript.Data>} object on
2656 L{Function<escript.Function>}
2657 @keyword A_reduced: value for coefficient C{A_reduced}
2658 @type A_reduced: any type that can be cast to a L{Data<escript.Data>}
2659 object on L{ReducedFunction<escript.ReducedFunction>}
2660 @keyword B: value for coefficient C{B}
2661 @type B: any type that can be cast to a L{Data<escript.Data>} object on
2662 L{Function<escript.Function>}
2663 @keyword B_reduced: value for coefficient C{B_reduced}
2664 @type B_reduced: any type that can be cast to a L{Data<escript.Data>}
2665 object on L{ReducedFunction<escript.ReducedFunction>}
2666 @keyword C: value for coefficient C{C}
2667 @type C: any type that can be cast to a L{Data<escript.Data>} object on
2668 L{Function<escript.Function>}
2669 @keyword C_reduced: value for coefficient C{C_reduced}
2670 @type C_reduced: any type that can be cast to a L{Data<escript.Data>}
2671 object on L{ReducedFunction<escript.ReducedFunction>}
2672 @keyword D: value for coefficient C{D}
2673 @type D: any type that can be cast to a L{Data<escript.Data>} object on
2674 L{Function<escript.Function>}
2675 @keyword D_reduced: value for coefficient C{D_reduced}
2676 @type D_reduced: any type that can be cast to a L{Data<escript.Data>}
2677 object on L{ReducedFunction<escript.ReducedFunction>}
2678 @keyword X: value for coefficient C{X}
2679 @type X: any type that can be cast to a L{Data<escript.Data>} object on
2680 L{Function<escript.Function>}
2681 @keyword X_reduced: value for coefficient C{X_reduced}
2682 @type X_reduced: any type that can be cast to a L{Data<escript.Data>}
2683 object on L{ReducedFunction<escript.ReducedFunction>}
2684 @keyword Y: value for coefficient C{Y}
2685 @type Y: any type that can be cast to a L{Data<escript.Data>} object on
2686 L{Function<escript.Function>}
2687 @keyword Y_reduced: value for coefficient C{Y_reduced}
2688 @type Y_reduced: any type that can be cast to a L{Data<escript.Data>}
2689 object on L{ReducedFunction<escript.Function>}
2690 @keyword d: value for coefficient C{d}
2691 @type d: any type that can be cast to a L{Data<escript.Data>} object on
2692 L{FunctionOnBoundary<escript.FunctionOnBoundary>}
2693 @keyword d_reduced: value for coefficient C{d_reduced}
2694 @type d_reduced: any type that can be cast to a L{Data<escript.Data>}
2695 object on L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}
2696 @keyword y: value for coefficient C{y}
2697 @type y: any type that can be cast to a L{Data<escript.Data>} object on
2698 L{FunctionOnBoundary<escript.FunctionOnBoundary>}
2699 @keyword d_contact: value for coefficient C{d_contact}
2700 @type d_contact: any type that can be cast to a L{Data<escript.Data>}
2701 object on L{FunctionOnContactOne<escript.FunctionOnContactOne>}
2702 or L{FunctionOnContactZero<escript.FunctionOnContactZero>}
2703 @keyword d_contact_reduced: value for coefficient C{d_contact_reduced}
2704 @type d_contact_reduced: any type that can be cast to a L{Data<escript.Data>}
2705 object on L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}
2706 or L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>}
2707 @keyword y_contact: value for coefficient C{y_contact}
2708 @type y_contact: any type that can be cast to a L{Data<escript.Data>}
2709 object on L{FunctionOnContactOne<escript.FunctionOnContactOne>}
2710 or L{FunctionOnContactZero<escript.FunctionOnContactZero>}
2711 @keyword y_contact_reduced: value for coefficient C{y_contact_reduced}
2712 @type y_contact_reduced: any type that can be cast to a L{Data<escript.Data>}
2713 object on L{ReducedFunctionOnContactOne<escript.FunctionOnContactOne>}
2714 or L{ReducedFunctionOnContactZero<escript.FunctionOnContactZero>}
2715 @keyword r: values prescribed to the solution at the locations of
2716 constraints
2717 @type r: 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 solution
2720 @keyword q: mask for location of constraints
2721 @type q: any type that can be cast to a L{Data<escript.Data>} object on
2722 L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2723 depending on whether reduced order is used for the
2724 representation of the equation
2725 @raise IllegalCoefficient: if an unknown coefficient keyword is used
2726 """
2727 super(LinearPDE,self).setValue(**coefficients)
2728 # check if the systrem is inhomogeneous:
2729 if len(coefficients)>0 and not self.isUsingLumping():
2730 q=self.getCoefficient("q")
2731 r=self.getCoefficient("r")
2732 if not q.isEmpty() and not r.isEmpty():
2733 if util.Lsup(q*r)>0.:
2734 self.trace("Inhomogeneous constraint detected.")
2735 self.invalidateSystem()
2736
2737
2738 def getResidual(self,u=None):
2739 """
2740 Returns the residual of u or the current solution if u is not present.
2741
2742 @param u: argument in the residual calculation. It must be representable
2743 in L{self.getFunctionSpaceForSolution()}. If u is not present
2744 or equals C{None} the current solution is used.
2745 @type u: L{Data<escript.Data>} or None
2746 @return: residual of u
2747 @rtype: L{Data<escript.Data>}
2748 """
2749 if u==None:
2750 return self.getOperator()*self.getSolution()-self.getRightHandSide()
2751 else:
2752 return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())-self.getRightHandSide()
2753
2754 def getFlux(self,u=None):
2755 """
2756 Returns the flux M{J} for a given M{u}.
2757
2758 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]}
2759
2760 or
2761
2762 M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[l]+(B[j]+B_reduced[j])u-X[j]-X_reduced[j]}
2763
2764 @param u: argument in the flux. If u is not present or equals L{None} the
2765 current solution is used.
2766 @type u: L{Data<escript.Data>} or None
2767 @return: flux
2768 @rtype: L{Data<escript.Data>}
2769 """
2770 if u==None: u=self.getSolution()
2771 return util.tensormult(self.getCoefficient("A"),util.grad(u,Funtion(self.getDomain))) \
2772 +util.matrixmult(self.getCoefficient("B"),u) \
2773 -util.self.getCoefficient("X") \
2774 +util.tensormult(self.getCoefficient("A_reduced"),util.grad(u,ReducedFuntion(self.getDomain))) \
2775 +util.matrixmult(self.getCoefficient("B_reduced"),u) \
2776 -util.self.getCoefficient("X_reduced")
2777
2778
2779 class Poisson(LinearPDE):
2780 """
2781 Class to define a Poisson equation problem. This is generally a
2782 L{LinearPDE} of the form
2783
2784 M{-grad(grad(u)[j])[j] = f}
2785
2786 with natural boundary conditions
2787
2788 M{n[j]*grad(u)[j] = 0 }
2789
2790 and constraints:
2791
2792 M{u=0} where M{q>0}
2793
2794 """
2795
2796 def __init__(self,domain,debug=False):
2797 """
2798 Initializes a new Poisson equation.
2799
2800 @param domain: domain of the PDE
2801 @type domain: L{Domain<escript.Domain>}
2802 @param debug: if True debug information is printed
2803
2804 """
2805 super(Poisson, self).__init__(domain,1,1,debug)
2806 self.introduceCoefficients(
2807 f=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2808 f_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE))
2809 self.setSymmetryOn()
2810
2811 def setValue(self,**coefficients):
2812 """
2813 Sets new values to coefficients.
2814
2815 @param coefficients: new values assigned to coefficients
2816 @keyword f: value for right hand side M{f}
2817 @type f: any type that can be cast to a L{Scalar<escript.Scalar>} object
2818 on L{Function<escript.Function>}
2819 @keyword q: mask for location of constraints
2820 @type q: any type that can be cast to a rank zero L{Data<escript.Data>}
2821 object on L{Solution<escript.Solution>} or
2822 L{ReducedSolution<escript.ReducedSolution>} depending on whether
2823 reduced order is used for the representation of the equation
2824 @raise IllegalCoefficient: if an unknown coefficient keyword is used
2825 """
2826 super(Poisson, self).setValue(**coefficients)
2827
2828
2829 def getCoefficient(self,name):
2830 """
2831 Returns the value of the coefficient C{name} of the general PDE.
2832
2833 @param name: name of the coefficient requested
2834 @type name: C{string}
2835 @return: the value of the coefficient C{name}
2836 @rtype: L{Data<escript.Data>}
2837 @raise IllegalCoefficient: invalid coefficient name
2838 @note: This method is called by the assembling routine to map the Poisson
2839 equation onto the general PDE.
2840 """
2841 if name == "A" :
2842 return escript.Data(util.kronecker(self.getDim()),escript.Function(self.getDomain()))
2843 elif name == "Y" :
2844 return self.getCoefficient("f")
2845 elif name == "Y_reduced" :
2846 return self.getCoefficient("f_reduced")
2847 else:
2848 return super(Poisson, self).getCoefficient(name)
2849
2850 class Helmholtz(LinearPDE):
2851 """
2852 Class to define a Helmholtz equation problem. This is generally a
2853 L{LinearPDE} of the form
2854
2855 M{S{omega}*u - grad(k*grad(u)[j])[j] = f}
2856
2857 with natural boundary conditions
2858
2859 M{k*n[j]*grad(u)[j] = g- S{alpha}u }
2860
2861 and constraints:
2862
2863 M{u=r} where M{q>0}
2864
2865 """
2866
2867 def __init__(self,domain,debug=False):
2868 """
2869 Initializes a new Helmholtz equation.
2870
2871 @param domain: domain of the PDE
2872 @type domain: L{Domain<escript.Domain>}
2873 @param debug: if True debug information is printed
2874
2875 """
2876 super(Helmholtz, self).__init__(domain,1,1,debug)
2877 self.introduceCoefficients(
2878 omega=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.OPERATOR),
2879 k=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.OPERATOR),
2880 f=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2881 f_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2882 alpha=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.OPERATOR),
2883 g=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2884 g_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE))
2885 self.setSymmetryOn()
2886
2887 def setValue(self,**coefficients):
2888 """
2889 Sets new values to coefficients.
2890
2891 @param coefficients: new values assigned to coefficients
2892 @keyword omega: value for coefficient M{S{omega}}
2893 @type omega: any type that can be cast to a L{Scalar<escript.Scalar>}
2894 object on L{Function<escript.Function>}
2895 @keyword k: value for coefficient M{k}
2896 @type k: any type that can be cast to a L{Scalar<escript.Scalar>} object
2897 on L{Function<escript.Function>}
2898 @keyword f: value for right hand side M{f}
2899 @type f: any type that can be cast to a L{Scalar<escript.Scalar>} object
2900 on L{Function<escript.Function>}
2901 @keyword alpha: value for right hand side M{S{alpha}}
2902 @type alpha: any type that can be cast to a L{Scalar<escript.Scalar>}
2903 object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}
2904 @keyword g: value for right hand side M{g}
2905 @type g: any type that can be cast to a L{Scalar<escript.Scalar>} object
2906 on L{FunctionOnBoundary<escript.FunctionOnBoundary>}
2907 @keyword r: prescribed values M{r} for the solution in constraints
2908 @type r: any type that can be cast to a L{Scalar<escript.Scalar>} object
2909 on L{Solution<escript.Solution>} or
2910 L{ReducedSolution<escript.ReducedSolution>} depending on whether
2911 reduced order is used for the representation of the equation
2912 @keyword q: mask for the location of constraints
2913 @type q: any type that can be cast to a L{Scalar<escript.Scalar>} object
2914 on L{Solution<escript.Solution>} or
2915 L{ReducedSolution<escript.ReducedSolution>} depending on whether
2916 reduced order is used for the representation of the equation
2917 @raise IllegalCoefficient: if an unknown coefficient keyword is used
2918 """
2919 super(Helmholtz, self).setValue(**coefficients)
2920
2921 def getCoefficient(self,name):
2922 """
2923 Returns the value of the coefficient C{name} of the general PDE.
2924
2925 @param name: name of the coefficient requested
2926 @type name: C{string}
2927 @return: the value of the coefficient C{name}
2928 @rtype: L{Data<escript.Data>}
2929 @raise IllegalCoefficient: invalid name
2930 """
2931 if name == "A" :
2932 if self.getCoefficient("k").isEmpty():
2933 return escript.Data(numpy.identity(self.getDim()),escript.Function(self.getDomain()))
2934 else:
2935 return escript.Data(numpy.identity(self.getDim()),escript.Function(self.getDomain()))*self.getCoefficient("k")
2936 elif name == "D" :
2937 return self.getCoefficient("omega")
2938 elif name == "Y" :
2939 return self.getCoefficient("f")
2940 elif name == "d" :
2941 return self.getCoefficient("alpha")
2942 elif name == "y" :
2943 return self.getCoefficient("g")
2944 elif name == "Y_reduced" :
2945 return self.getCoefficient("f_reduced")
2946 elif name == "y_reduced" :
2947 return self.getCoefficient("g_reduced")
2948 else:
2949 return super(Helmholtz, self).getCoefficient(name)
2950
2951 class LameEquation(LinearPDE):
2952 """
2953 Class to define a Lame equation problem. This problem is defined as:
2954
2955 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] }
2956
2957 with natural boundary conditions:
2958
2959 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] }
2960
2961 and constraints:
2962
2963 M{u[i]=r[i]} where M{q[i]>0}
2964
2965 """
2966
2967 def __init__(self,domain,debug=False):
2968 """
2969 Initializes a new Lame equation.
2970
2971 @param domain: domain of the PDE
2972 @type domain: L{Domain<escript.Domain>}
2973 @param debug: if True debug information is printed
2974
2975 """
2976 super(LameEquation, self).__init__(domain,\
2977 domain.getDim(),domain.getDim(),debug)
2978 self.introduceCoefficients(lame_lambda=PDECoef(PDECoef.INTERIOR,(),PDECoef.OPERATOR),
2979 lame_mu=PDECoef(PDECoef.INTERIOR,(),PDECoef.OPERATOR),
2980 F=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2981 sigma=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
2982 f=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE))
2983 self.setSymmetryOn()
2984
2985 def setValues(self,**coefficients):
2986 """
2987 Sets new values to coefficients.
2988
2989 @param coefficients: new values assigned to coefficients
2990 @keyword lame_mu: value for coefficient M{S{mu}}
2991 @type lame_mu: any type that can be cast to a L{Scalar<escript.Scalar>}
2992 object on L{Function<escript.Function>}
2993 @keyword lame_lambda: value for coefficient M{S{lambda}}
2994 @type lame_lambda: any type that can be cast to a L{Scalar<escript.Scalar>}
2995 object on L{Function<escript.Function>}
2996 @keyword F: value for internal force M{F}
2997 @type F: any type that can be cast to a L{Vector<escript.Vector>} object
2998 on L{Function<escript.Function>}
2999 @keyword sigma: value for initial stress M{S{sigma}}
3000 @type sigma: any type that can be cast to a L{Tensor<escript.Tensor>}
3001 object on L{Function<escript.Function>}
3002 @keyword f: value for external force M{f}
3003 @type f: any type that can be cast to a L{Vector<escript.Vector>} object
3004 on L{FunctionOnBoundary<escript.FunctionOnBoundary>}
3005 @keyword r: prescribed values M{r} for the solution in constraints
3006 @type r: any type that can be cast to a L{Vector<escript.Vector>} object
3007 on L{Solution<escript.Solution>} or
3008 L{ReducedSolution<escript.ReducedSolution>} depending on whether
3009 reduced order is used for the representation of the equation
3010 @keyword q: mask for the location of constraints
3011 @type q: any type that can be cast to a L{Vector<escript.Vector>} object
3012 on L{Solution<escript.Solution>} or
3013 L{ReducedSolution<escript.ReducedSolution>} depending on whether
3014 reduced order is used for the representation of the equation
3015 @raise IllegalCoefficient: if an unknown coefficient keyword is used
3016 """
3017 super(LameEquation, self).setValues(**coefficients)
3018
3019 def getCoefficient(self,name):
3020 """
3021 Returns the value of the coefficient C{name} of the general PDE.
3022
3023 @param name: name of the coefficient requested
3024 @type name: C{string}
3025 @return: the value of the coefficient C{name}
3026 @rtype: L{Data<escript.Data>}
3027 @raise IllegalCoefficient: invalid coefficient name
3028 """
3029 out =self.createCoefficient("A")
3030 if name == "A" :
3031 if self.getCoefficient("lame_lambda").isEmpty():
3032 if self.getCoefficient("lame_mu").isEmpty():
3033 pass
3034 else:
3035 for i in range(self.getDim()):
3036 for j in range(self.getDim()):
3037 out[i,j,j,i] += self.getCoefficient("lame_mu")
3038 out[i,j,i,j] += self.getCoefficient("lame_mu")
3039 else:
3040 if self.getCoefficient("lame_mu").isEmpty():
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 else:
3045 for i in range(self.getDim()):
3046 for j in range(self.getDim()):
3047 out[i,i,j,j] += self.getCoefficient("lame_lambda")
3048 out[i,j,j,i] += self.getCoefficient("lame_mu")
3049 out[i,j,i,j] += self.getCoefficient("lame_mu")
3050 return out
3051 elif name == "X" :
3052 return self.getCoefficient("sigma")
3053 elif name == "Y" :
3054 return self.getCoefficient("F")
3055 elif name == "y" :
3056 return self.getCoefficient("f")
3057 else:
3058 return super(LameEquation, self).getCoefficient(name)
3059
3060 def LinearSinglePDE(domain,debug=False):
3061 """
3062 Defines a single linear PDE.
3063
3064 @param domain: domain of the PDE
3065 @type domain: L{Domain<escript.Domain>}
3066 @param debug: if True debug information is printed
3067 @rtype: L{LinearPDE}
3068 """
3069 return LinearPDE(domain,numEquations=1,numSolutions=1,debug=debug)
3070
3071 def LinearPDESystem(domain,debug=False):
3072 """
3073 Defines a system of linear PDEs.
3074
3075 @param domain: domain of the PDEs
3076 @type domain: L{Domain<escript.Domain>}
3077 @param debug: if True debug information is printed
3078 @rtype: L{LinearPDE}
3079 """
3080 return LinearPDE(domain,numEquations=domain.getDim(),numSolutions=domain.getDim(),debug=debug)
3081
3082
3083 class TransportPDE(LinearProblem):
3084 """
3085 This class is used to define a transport problem given by a general linear,
3086 time dependent, second order PDE for an unknown, non-negative,
3087 time-dependent function M{u} on a given domain defined through a
3088 L{Domain<escript.Domain>} object.
3089
3090 For a single equation with a solution with a single component the transport
3091 problem is defined in the following form:
3092
3093 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)}
3094
3095 where M{u_t} denotes the time derivative of M{u} and M{grad(F)} denotes the
3096 spatial derivative of M{F}. Einstein's summation convention, ie. summation
3097 over indexes appearing twice in a term of a sum performed, is used.
3098 The coefficients M{M}, M{A}, M{B}, M{C}, M{D}, M{X} and M{Y} have to be
3099 specified through L{Data<escript.Data>} objects in L{Function<escript.Function>}
3100 and the coefficients M{M_reduced}, M{A_reduced}, M{B_reduced}, M{C_reduced},
3101 M{D_reduced}, M{X_reduced} and M{Y_reduced} have to be specified through
3102 L{Data<escript.Data>} objects in L{ReducedFunction<escript.ReducedFunction>}.
3103 It is also allowed to use objects that can be converted into such
3104 L{Data<escript.Data>} objects. M{M} and M{M_reduced} are scalar, M{A} and
3105 M{A_reduced} are rank two, M{B}, M{C}, M{X}, M{B_reduced}, M{C_reduced} and
3106 M{X_reduced} are rank one and M{D}, M{D_reduced}, M{Y} and M{Y_reduced} are
3107 scalar.
3108
3109 The following natural boundary conditions are considered:
3110
3111 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}
3112
3113 where M{n} is the outer normal field. Notice that the coefficients M{A},
3114 M{A_reduced}, M{B}, M{B_reduced}, M{X} and M{X_reduced} are defined in the
3115 transport problem. The coefficients M{m}, M{d} and M{y} are each a scalar in
3116 L{FunctionOnBoundary<escript.FunctionOnBoundary>} and the coefficients
3117 M{m_reduced}, M{d_reduced} and M{y_reduced} are each a scalar in
3118 L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
3119
3120 Constraints for the solution prescribing the value of the solution at
3121 certain locations in the domain have the form
3122
3123 M{u_t=r} where M{q>0}
3124
3125 M{r} and M{q} are each scalar where M{q} is the characteristic function
3126 defining where the constraint is applied. The constraints override any other
3127 condition set by the transport problem or the boundary condition.
3128
3129 The transport problem is symmetrical if
3130
3131 M{A[i,j]=A[j,i]} and M{B[j]=C[j]} and M{A_reduced[i,j]=A_reduced[j,i]} and
3132 M{B_reduced[j]=C_reduced[j]}
3133
3134 For a system and a solution with several components the transport problem
3135 has the form
3136
3137 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] }
3138
3139 M{A} and M{A_reduced} are of rank four, M{B}, M{B_reduced}, M{C} and
3140 M{C_reduced} are each of rank three, M{M}, M{M_reduced}, M{D}, M{D_reduced},
3141 M{X_reduced} and M{X} are each of rank two and M{Y} and M{Y_reduced} are of
3142 rank one. The natural boundary conditions take the form:
3143
3144 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}
3145
3146 The coefficient M{d} and M{m} are of rank two and M{y} is of rank one with
3147 all in L{FunctionOnBoundary<escript.FunctionOnBoundary>}. The coefficients
3148 M{d_reduced} and M{m_reduced} are of rank two and M{y_reduced} is of rank
3149 one all in L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
3150
3151 Constraints take the form
3152
3153 M{u[i]_t=r[i]} where M{q[i]>0}
3154
3155 M{r} and M{q} are each rank one. Notice that at some locations not
3156 necessarily all components must have a constraint.
3157
3158 The transport problem is symmetrical if
3159
3160 - M{M[i,k]=M[i,k]}
3161 - M{M_reduced[i,k]=M_reduced[i,k]}
3162 - M{A[i,j,k,l]=A[k,l,i,j]}
3163 - M{A_reduced[i,j,k,l]=A_reduced[k,l,i,j]}
3164 - M{B[i,j,k]=C[k,i,j]}
3165 - M{B_reduced[i,j,k]=C_reduced[k,i,j]}
3166 - M{D[i,k]=D[i,k]}
3167 - M{D_reduced[i,k]=D_reduced[i,k]}
3168 - M{m[i,k]=m[k,i]}
3169 - M{m_reduced[i,k]=m_reduced[k,i]}
3170 - M{d[i,k]=d[k,i]}
3171 - M{d_reduced[i,k]=d_reduced[k,i]}
3172
3173 L{TransportPDE} also supports solution discontinuities over a contact region
3174 in the domain. To specify the conditions across the discontinuity we are
3175 using the generalised flux M{J} which, in the case of a system of PDEs and
3176 several components of the solution, is defined as
3177
3178 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]}
3179
3180 For the case of single solution component and single PDE M{J} is defined as
3181
3182 M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[j]+(B[i]+B_reduced[i])*u+X[i]+X_reduced[i]}
3183
3184 In the context of discontinuities M{n} denotes the normal on the
3185 discontinuity pointing from side 0 towards side 1 calculated from
3186 L{getNormal<escript.FunctionSpace.getNormal>} of L{FunctionOnContactZero<escript.FunctionOnContactZero>}.
3187 For a system of transport problems the contact condition takes the form
3188
3189 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]}
3190
3191 where M{J0} and M{J1} are the fluxes on side 0 and side 1 of the
3192 discontinuity, respectively. M{jump(u)}, which is the difference of the
3193 solution at side 1 and at side 0, denotes the jump of M{u} across
3194 discontinuity along the normal calculated by L{jump<util.jump>}.
3195 The coefficient M{d_contact} is of rank two and M{y_contact} is of rank one
3196 both in L{FunctionOnContactZero<escript.FunctionOnContactZero>} or L{FunctionOnContactOne<escript.FunctionOnContactOne>}.
3197 The coefficient M{d_contact_reduced} is of rank two and M{y_contact_reduced}
3198 is of rank one both in L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}.
3199 In case of a single PDE and a single component solution the contact
3200 condition takes the form
3201
3202 M{n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)}
3203
3204 In this case the coefficient M{d_contact} and M{y_contact} are each scalar
3205 both in L{FunctionOnContactZero<escript.FunctionOnContactZero>} or
3206 L{FunctionOnContactOne<escript.FunctionOnContactOne>} and the coefficient
3207 M{d_contact_reduced} and M{y_contact_reduced} are each scalar both in
3208 L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or
3209 L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}.
3210
3211 Typical usage::
3212
3213 p = TransportPDE(dom)
3214 p.setValue(M=1., C=[-1.,0.])
3215 p.setInitialSolution(u=exp(-length(dom.getX()-[0.1,0.1])**2)
3216 t = 0
3217 dt = 0.1
3218 while (t < 1.):
3219 u = p.solve(dt)
3220
3221 """
3222 def __init__(self,domain,numEquations=None,numSolutions=None, useBackwardEuler=False, debug=False):
3223 """
3224 Initializes a transport problem.
3225
3226 @param domain: domain of the PDE
3227 @type domain: L{Domain<escript.Domain>}
3228 @param numEquations: number of equations. If C{None} the number of
3229 equations is extracted from the coefficients.
3230 @param numSolutions: number of solution components. If C{None} the number
3231 of solution components is extracted from the
3232 coefficients.
3233 @param debug: if True debug information is printed
3234 @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.
3235 @type useBackwardEuler: C{bool}
3236 """
3237 if useBackwardEuler:
3238 self.__useBackwardEuler=True
3239 else:
3240 self.__useBackwardEuler=False
3241 super(TransportPDE, self).__init__(domain,numEquations,numSolutions,debug)
3242
3243 self.setConstraintWeightingFactor()
3244 #
3245 # the coefficients of the transport problem
3246 #
3247 self.introduceCoefficients(
3248 M=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
3249 A=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
3250 B=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
3251 C=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
3252 D=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
3253 X=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
3254 Y=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
3255 m=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
3256 d=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
3257 y=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
3258 d_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
3259 y_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
3260 M_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
3261 A_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
3262 B_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
3263 C_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
3264 D_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
3265 X_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
3266 Y_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
3267 m_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
3268 d_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
3269 y_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
3270 d_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
3271 y_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
3272 r=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.RIGHTHANDSIDE),
3273 q=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.BOTH) )
3274
3275 def __str__(self):
3276 """
3277 Returns the string representation of the transport problem.
3278
3279 @return: a simple representation of the transport problem
3280 @rtype: C{str}
3281 """
3282 return "<TransportPDE %d>"%id(self)
3283
3284 def useBackwardEuler(self):
3285 """
3286 Returns true if backward Euler is used. Otherwise false is returned.
3287 @rtype: bool
3288 """
3289 return self.__useBackwardEuler
3290
3291
3292 def checkSymmetry(self,verbose=True):
3293 """
3294 Tests the transport problem for symmetry.
3295
3296 @param verbose: if set to True or not present a report on coefficients
3297 which break the symmetry is printed.
3298 @type verbose: C{bool}
3299 @return: True if the PDE is symmetric
3300 @rtype: C{bool}
3301 @note: This is a very expensive operation. It should be used for
3302 degugging only! The symmetry flag is not altered.
3303 """
3304 out=True
3305 out=out and self.checkSymmetricTensor("M", verbose)
3306 out=out and self.checkSymmetricTensor("M_reduced", verbose)
3307 out=out and self.checkSymmetricTensor("A", verbose)
3308 out=out and self.checkSymmetricTensor("A_reduced", verbose)
3309 out=out and self.checkReciprocalSymmetry("B","C", verbose)
3310 out=out and self.checkReciprocalSymmetry("B_reduced","C_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("m", verbose)
3314 out=out and self.checkSymmetricTensor("m_reduced", verbose)
3315 out=out and self.checkSymmetricTensor("d", verbose)
3316 out=out and self.checkSymmetricTensor("d_reduced", verbose)
3317 out=out and self.checkSymmetricTensor("d_contact", verbose)
3318 out=out and self.checkSymmetricTensor("d_contact_reduced", verbose)
3319 return out
3320
3321 def setValue(self,**coefficients):
3322 """
3323 Sets new values to coefficients.
3324
3325 @param coefficients: new values assigned to coefficients
3326 @keyword M: value for coefficient C{M}
3327 @type M: any type that can be cast to a L{Data<escript.Data>} object on
3328 L{Function<escript.Function>}
3329 @keyword M_reduced: value for coefficient C{M_reduced}
3330 @type M_reduced: any type that can be cast to a L{Data<escript.Data>}
3331 object on L{Function<escript.ReducedFunction>}
3332 @keyword A: value for coefficient C{A}
3333 @type A: any type that can be cast to a L{Data<escript.Data>} object on
3334 L{Function<escript.Function>}
3335 @keyword A_reduced: value for coefficient C{A_reduced}
3336 @type A_reduced: any type that can be cast to a L{Data<escript.Data>}
3337 object on L{ReducedFunction<escript.ReducedFunction>}
3338 @keyword B: value for coefficient C{B}
3339 @type B: any type that can be cast to a L{Data<escript.Data>} object on
3340 L{Function<escript.Function>}
3341 @keyword B_reduced: value for coefficient C{B_reduced}
3342 @type B_reduced: any type that can be cast to a L{Data<escript.Data>}
3343 object on L{ReducedFunction<escript.ReducedFunction>}
3344 @keyword C: value for coefficient C{C}
3345 @type C: any type that can be cast to a L{Data<escript.Data>} object on
3346 L{Function<escript.Function>}
3347 @keyword C_reduced: value for coefficient C{C_reduced}
3348 @type C_reduced: any type that can be cast to a L{Data<escript.Data>}
3349 object on L{ReducedFunction<escript.ReducedFunction>}
3350 @keyword D: value for coefficient C{D}
3351 @type D: any type that can be cast to a L{Data<escript.Data>} object on
3352 L{Function<escript.Function>}
3353 @keyword D_reduced: value for coefficient C{D_reduced}
3354 @type D_reduced: any type that can be cast to a L{Data<escript.Data>}
3355 object on L{ReducedFunction<escript.ReducedFunction>}
3356 @keyword X: value for coefficient C{X}
3357 @type X: any type that can be cast to a L{Data<escript.Data>} object on
3358 L{Function<escript.Function>}
3359 @keyword X_reduced: value for coefficient C{X_reduced}
3360 @type X_reduced: any type that can be cast to a L{Data<escript.Data>}
3361 object on L{ReducedFunction<escript.ReducedFunction>}
3362 @keyword Y: value for coefficient C{Y}
3363 @type Y: any type that can be cast to a L{Data<escript.Data>} object on
3364 L{Function<escript.Function>}
3365 @keyword Y_reduced: value for coefficient C{Y_reduced}
3366 @type Y_reduced: any type that can be cast to a L{Data<escript.Data>}
3367 object on L{ReducedFunction<escript.Function>}
3368 @keyword m: value for coefficient C{m}
3369 @type m: any type that can be cast to a L{Data<escript.Data>} object on
3370 L{FunctionOnBoundary<escript.FunctionOnBoundary>}
3371 @keyword m_reduced: value for coefficient C{m_reduced}
3372 @type m_reduced: any type that can be cast to a L{Data<escript.Data>}
3373 object on L{FunctionOnBoundary<escript.ReducedFunctionOnBoundary>}
3374 @keyword d: value for coefficient C{d}
3375 @type d: any type that can be cast to a L{Data<escript.Data>} object on
3376 L{FunctionOnBoundary<escript.FunctionOnBoundary>}
3377 @keyword d_reduced: value for coefficient C{d_reduced}
3378 @type d_reduced: any type that can be cast to a L{Data<escript.Data>}
3379 object on L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}
3380 @keyword y: value for coefficient C{y}
3381 @type y: any type that can be cast to a L{Data<escript.Data>} object on
3382 L{FunctionOnBoundary<escript.FunctionOnBoundary>}
3383 @keyword d_contact: value for coefficient C{d_contact}
3384 @type d_contact: any type that can be cast to a L{Data<escript.Data>}
3385 object on L{FunctionOnContactOne<escript.FunctionOnContactOne>} or L{FunctionOnContactZero<escript.FunctionOnContactZero>}
3386 @keyword d_contact_reduced: value for coefficient C{d_contact_reduced}
3387 @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>}
3388 @keyword y_contact: value for coefficient C{y_contact}
3389 @type y_contact: any type that can be cast to a L{Data<escript.Data>}
3390 object on L{FunctionOnContactOne<escript.FunctionOnContactOne>} or L{FunctionOnContactZero<escript.FunctionOnContactZero>}
3391 @keyword y_contact_reduced: value for coefficient C{y_contact_reduced}
3392 @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>}
3393 @keyword r: values prescribed to the solution at the locations of constraints
3394 @type r: any type that can be cast to a L{Data<escript.Data>} object on
3395 L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
3396 depending on whether reduced order is used for the solution
3397 @keyword q: mask for the location of constraints
3398 @type q: any type that can be cast to a L{Data<escript.Data>} object on
3399 L{Solution<escript.Solution>} or
3400 L{ReducedSolution<escript.ReducedSolution>} depending on whether
3401 reduced order is used for the representation of the equation
3402 @raise IllegalCoefficient: if an unknown coefficient keyword is used
3403 """
3404 super(TransportPDE,self).setValue(**coefficients)
3405
3406 def createOperator(self):
3407 """
3408 Returns an instance of a new transport operator.
3409 """
3410 if self.useBackwardEuler():
3411 theta=1.
3412 else:
3413 theta=