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