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