/[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 3094 - (show annotations)
Fri Aug 13 08:38:06 2010 UTC (8 years, 11 months ago) by gross
File MIME type: text/x-python
File size: 153310 byte(s)
The MPI and sequational GAUSS_SEIDEL have been merged.
The couring and main diagonal pointer is now manged by the patternm which means that they are calculated once only even if the preconditioner is deleted.



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