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