/[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 3283 - (show annotations)
Mon Oct 18 22:39:28 2010 UTC (10 years, 4 months ago) by gross
File MIME type: text/x-python
File size: 156835 byte(s)
AMG reengineered, BUG is Smoother fixed.


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