/[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 3259 - (show annotations)
Mon Oct 11 01:48:14 2010 UTC (8 years, 11 months ago) by jfenwick
File MIME type: text/x-python
File size: 154016 byte(s)
Merging dudley and scons updates from branches

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