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