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