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