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