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