/[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 1841 - (show annotations)
Fri Oct 3 03:57:52 2008 UTC (11 years, 1 month ago) by gross
File MIME type: text/x-python
File size: 115272 byte(s)
modification on LinearPDE class and a first version of Transport class
1
2 ########################################################
3 #
4 # Copyright (c) 2003-2008 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-2008 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__="http://www.uq.edu.au/esscc/escript-finley"
21
22 """
23 The module provides an interface to define and solve linear partial
24 differential equations (PDEs) and Transport problems within L{escript}. L{linearPDEs} does not provide any
25 solver capabilities in itself but hands the PDE over to
26 the PDE solver library defined through the L{Domain<escript.Domain>} of the PDE.
27 The general interface is provided through the L{LinearPDE} class.
28 L{TransportProblem} provides an interface to initial value problems dominated by its advective terms.
29
30 @var __author__: name of author
31 @var __copyright__: copyrights
32 @var __license__: licence agreement
33 @var __url__: url entry point on documentation
34 @var __version__: version
35 @var __date__: date of the version
36 """
37
38 import math
39 import escript
40 import util
41 import numarray
42
43 __author__="Lutz Gross, l.gross@uq.edu.au"
44
45
46 class IllegalCoefficient(ValueError):
47 """
48 raised if an illegal coefficient of the general ar particular PDE is requested.
49 """
50 pass
51
52 class IllegalCoefficientValue(ValueError):
53 """
54 raised if an incorrect value for a coefficient is used.
55 """
56 pass
57
58 class IllegalCoefficientFunctionSpace(ValueError):
59 """
60 raised if an incorrect function space for a coefficient is used.
61 """
62
63 class UndefinedPDEError(ValueError):
64 """
65 raised if a PDE is not fully defined yet.
66 """
67 pass
68
69 class PDECoef(object):
70 """
71 A class for describing a PDE coefficient
72
73 @cvar INTERIOR: indicator that coefficient is defined on the interior of the PDE domain
74 @cvar BOUNDARY: indicator that coefficient is defined on the boundary of the PDE domain
75 @cvar CONTACT: indicator that coefficient is defined on the contact region within the PDE domain
76 @cvar INTERIOR_REDUCED: indicator that coefficient is defined on the interior of the PDE domain using a reduced integration order
77 @cvar BOUNDARY_REDUCED: indicator that coefficient is defined on the boundary of the PDE domain using a reduced integration order
78 @cvar CONTACT_REDUCED: indicator that coefficient is defined on the contact region within the PDE domain using a reduced integration order
79 @cvar SOLUTION: indicator that coefficient is defined trough a solution of the PDE
80 @cvar REDUCED: indicator that coefficient is defined trough a reduced solution of the PDE
81 @cvar BY_EQUATION: indicator that the dimension of the coefficient shape is defined by the number PDE equations
82 @cvar BY_SOLUTION: indicator that the dimension of the coefficient shape is defined by the number PDE solutions
83 @cvar BY_DIM: indicator that the dimension of the coefficient shape is defined by the spatial dimension
84 @cvar OPERATOR: indicator that the the coefficient alters the operator of the PDE
85 @cvar RIGHTHANDSIDE: indicator that the the coefficient alters the right hand side of the PDE
86 @cvar BOTH: indicator that the the coefficient alters the operator as well as the right hand side of the PDE
87
88 """
89 INTERIOR=0
90 BOUNDARY=1
91 CONTACT=2
92 SOLUTION=3
93 REDUCED=4
94 BY_EQUATION=5
95 BY_SOLUTION=6
96 BY_DIM=7
97 OPERATOR=10
98 RIGHTHANDSIDE=11
99 BOTH=12
100 INTERIOR_REDUCED=13
101 BOUNDARY_REDUCED=14
102 CONTACT_REDUCED=15
103
104 def __init__(self, where, pattern, altering):
105 """
106 Initialise a PDE Coefficient type
107
108 @param where: describes where the coefficient lives
109 @type where: one of L{INTERIOR}, L{BOUNDARY}, L{CONTACT}, L{SOLUTION}, L{REDUCED},
110 L{INTERIOR_REDUCED}, L{BOUNDARY_REDUCED}, L{CONTACT_REDUCED}.
111 @param pattern: describes the shape of the coefficient and how the shape is build for a given
112 spatial dimension and numbers of equation and solution in then PDE. For instance,
113 (L{BY_EQUATION},L{BY_SOLUTION},L{BY_DIM}) descrbes a rank 3 coefficient which
114 is instanciated as shape (3,2,2) in case of a three equations and two solution components
115 on a 2-dimensional domain. In the case of single equation and a single solution component
116 the shape compoments marked by L{BY_EQUATION} or L{BY_SOLUTION} are dropped. In this case
117 the example would be read as (2,).
118 @type pattern: C{tuple} of L{BY_EQUATION}, L{BY_SOLUTION}, L{BY_DIM}
119 @param altering: indicates what part of the PDE is altered if the coefficiennt is altered
120 @type altering: one of L{OPERATOR}, L{RIGHTHANDSIDE}, L{BOTH}
121 @param reduced: indicates if reduced
122 @type reduced: C{bool}
123 """
124 super(PDECoef, self).__init__()
125 self.what=where
126 self.pattern=pattern
127 self.altering=altering
128 self.resetValue()
129
130 def resetValue(self):
131 """
132 resets coefficient value to default
133 """
134 self.value=escript.Data()
135
136 def getFunctionSpace(self,domain,reducedEquationOrder=False,reducedSolutionOrder=False):
137 """
138 defines the L{FunctionSpace<escript.FunctionSpace>} of the coefficient
139
140 @param domain: domain on which the PDE uses the coefficient
141 @type domain: L{Domain<escript.Domain>}
142 @param reducedEquationOrder: True to indicate that reduced order is used to represent the equation
143 @type reducedEquationOrder: C{bool}
144 @param reducedSolutionOrder: True to indicate that reduced order is used to represent the solution
145 @type reducedSolutionOrder: C{bool}
146 @return: L{FunctionSpace<escript.FunctionSpace>} of the coefficient
147 @rtype: L{FunctionSpace<escript.FunctionSpace>}
148 """
149 if self.what==self.INTERIOR:
150 return escript.Function(domain)
151 elif self.what==self.INTERIOR_REDUCED:
152 return escript.ReducedFunction(domain)
153 elif self.what==self.BOUNDARY:
154 return escript.FunctionOnBoundary(domain)
155 elif self.what==self.BOUNDARY_REDUCED:
156 return escript.ReducedFunctionOnBoundary(domain)
157 elif self.what==self.CONTACT:
158 return escript.FunctionOnContactZero(domain)
159 elif self.what==self.CONTACT_REDUCED:
160 return escript.ReducedFunctionOnContactZero(domain)
161 elif self.what==self.SOLUTION:
162 if reducedEquationOrder and reducedSolutionOrder:
163 return escript.ReducedSolution(domain)
164 else:
165 return escript.Solution(domain)
166 elif self.what==self.REDUCED:
167 return escript.ReducedSolution(domain)
168
169 def getValue(self):
170 """
171 returns the value of the coefficient
172
173 @return: value of the coefficient
174 @rtype: L{Data<escript.Data>}
175 """
176 return self.value
177
178 def setValue(self,domain,numEquations=1,numSolutions=1,reducedEquationOrder=False,reducedSolutionOrder=False,newValue=None):
179 """
180 set the value of the coefficient to a new value
181
182 @param domain: domain on which the PDE uses the coefficient
183 @type domain: L{Domain<escript.Domain>}
184 @param numEquations: number of equations of the PDE
185 @type numEquations: C{int}
186 @param numSolutions: number of components of the PDE solution
187 @type numSolutions: C{int}
188 @param reducedEquationOrder: True to indicate that reduced order is used to represent the equation
189 @type reducedEquationOrder: C{bool}
190 @param reducedSolutionOrder: True to indicate that reduced order is used to represent the solution
191 @type reducedSolutionOrder: C{bool}
192 @param newValue: number of components of the PDE solution
193 @type newValue: any object that can be converted into a L{Data<escript.Data>} object with the appropriate shape and L{FunctionSpace<escript.FunctionSpace>}
194 @raise IllegalCoefficientValue: if the shape of the assigned value does not match the shape of the coefficient
195 @raise IllegalCoefficientFunctionSpace: if unable to interploate value to appropriate function space
196 """
197 if newValue==None:
198 newValue=escript.Data()
199 elif isinstance(newValue,escript.Data):
200 if not newValue.isEmpty():
201 if not newValue.getFunctionSpace() == self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder):
202 try:
203 newValue=escript.Data(newValue,self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder))
204 except:
205 raise IllegalCoefficientFunctionSpace,"Unable to interpolate coefficient to function space %s"%self.getFunctionSpace(domain)
206 else:
207 newValue=escript.Data(newValue,self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder))
208 if not newValue.isEmpty():
209 if not self.getShape(domain,numEquations,numSolutions)==newValue.getShape():
210 raise IllegalCoefficientValue,"Expected shape of coefficient is %s but actual shape is %s."%(self.getShape(domain,numEquations,numSolutions),newValue.getShape())
211 self.value=newValue
212
213 def isAlteringOperator(self):
214 """
215 checks if the coefficient alters the operator of the PDE
216
217 @return: True if the operator of the PDE is changed when the coefficient is changed
218 @rtype: C{bool}
219 """
220 if self.altering==self.OPERATOR or self.altering==self.BOTH:
221 return not None
222 else:
223 return None
224
225 def isAlteringRightHandSide(self):
226 """
227 checks if the coefficeint alters the right hand side of the PDE
228
229 @rtype: C{bool}
230 @return: True if the right hand side of the PDE is changed when the coefficient is changed
231 """
232 if self.altering==self.RIGHTHANDSIDE or self.altering==self.BOTH:
233 return not None
234 else:
235 return None
236
237 def estimateNumEquationsAndNumSolutions(self,domain,shape=()):
238 """
239 tries to estimate the number of equations and number of solutions if the coefficient has the given shape
240
241 @param domain: domain on which the PDE uses the coefficient
242 @type domain: L{Domain<escript.Domain>}
243 @param shape: suggested shape of the coefficient
244 @type shape: C{tuple} of C{int} values
245 @return: the number of equations and number of solutions of the PDE is the coefficient has shape s.
246 If no appropriate numbers could be identified, C{None} is returned
247 @rtype: C{tuple} of two C{int} values or C{None}
248 """
249 dim=domain.getDim()
250 if len(shape)>0:
251 num=max(shape)+1
252 else:
253 num=1
254 search=[]
255 if self.definesNumEquation() and self.definesNumSolutions():
256 for u in range(num):
257 for e in range(num):
258 search.append((e,u))
259 search.sort(self.__CompTuple2)
260 for item in search:
261 s=self.getShape(domain,item[0],item[1])
262 if len(s)==0 and len(shape)==0:
263 return (1,1)
264 else:
265 if s==shape: return item
266 elif self.definesNumEquation():
267 for e in range(num,0,-1):
268 s=self.getShape(domain,e,0)
269 if len(s)==0 and len(shape)==0:
270 return (1,None)
271 else:
272 if s==shape: return (e,None)
273
274 elif self.definesNumSolutions():
275 for u in range(num,0,-1):
276 s=self.getShape(domain,0,u)
277 if len(s)==0 and len(shape)==0:
278 return (None,1)
279 else:
280 if s==shape: return (None,u)
281 return None
282 def definesNumSolutions(self):
283 """
284 checks if the coefficient allows to estimate the number of solution components
285
286 @return: True if the coefficient allows an estimate of the number of solution components
287 @rtype: C{bool}
288 """
289 for i in self.pattern:
290 if i==self.BY_SOLUTION: return True
291 return False
292
293 def definesNumEquation(self):
294 """
295 checks if the coefficient allows to estimate the number of equations
296
297 @return: True if the coefficient allows an estimate of the number of equations
298 @rtype: C{bool}
299 """
300 for i in self.pattern:
301 if i==self.BY_EQUATION: return True
302 return False
303
304 def __CompTuple2(self,t1,t2):
305 """
306 Compare two tuples of possible number of equations and number of solutions
307
308 @param t1: The first tuple
309 @param t2: The second tuple
310
311 """
312
313 dif=t1[0]+t1[1]-(t2[0]+t2[1])
314 if dif<0: return 1
315 elif dif>0: return -1
316 else: return 0
317
318 def getShape(self,domain,numEquations=1,numSolutions=1):
319 """
320 builds the required shape of the coefficient
321
322 @param domain: domain on which the PDE uses the coefficient
323 @type domain: L{Domain<escript.Domain>}
324 @param numEquations: number of equations of the PDE
325 @type numEquations: C{int}
326 @param numSolutions: number of components of the PDE solution
327 @type numSolutions: C{int}
328 @return: shape of the coefficient
329 @rtype: C{tuple} of C{int} values
330 """
331 dim=domain.getDim()
332 s=()
333 for i in self.pattern:
334 if i==self.BY_EQUATION:
335 if numEquations>1: s=s+(numEquations,)
336 elif i==self.BY_SOLUTION:
337 if numSolutions>1: s=s+(numSolutions,)
338 else:
339 s=s+(dim,)
340 return s
341
342 class LinearProblem(object):
343 """
344 This is the base class is to define a general linear PDE-type problem for an u
345 for an unknown function M{u} on a given domain defined through a L{Domain<escript.Domain>} object.
346 The problem can be given as a single equations or as a system of equations.
347
348 The class assumes that some sort of assembling process is required to form a problem of the form
349
350 M{L u=f}
351
352 where M{L} is an operator and M{f} is the right hand side. This operator problem will be solved to get the
353 unknown M{u}.
354
355 @cvar DEFAULT: The default method used to solve the system of linear equations
356 @cvar DIRECT: The direct solver based on LDU factorization
357 @cvar CHOLEVSKY: The direct solver based on LDLt factorization (can only be applied for symmetric PDEs)
358 @cvar PCG: The preconditioned conjugate gradient method (can only be applied for symmetric PDEs)
359 @cvar CR: The conjugate residual method
360 @cvar CGS: The conjugate gardient square method
361 @cvar BICGSTAB: The stabilized BiConjugate Gradient method.
362 @cvar TFQMR: Transport Free Quasi Minimal Residual method.
363 @cvar MINRES: Minimum residual method.
364 @cvar SSOR: The symmetric overrealaxtion method
365 @cvar ILU0: The incomplete LU factorization preconditioner with no fill in
366 @cvar ILUT: The incomplete LU factorization preconditioner with will in
367 @cvar JACOBI: The Jacobi preconditioner
368 @cvar GMRES: The Gram-Schmidt minimum residual method
369 @cvar PRES20: Special GMRES with restart after 20 steps and truncation after 5 residuals
370 @cvar LUMPING: Matrix lumping.
371 @cvar NO_REORDERING: No matrix reordering allowed
372 @cvar MINIMUM_FILL_IN: Reorder matrix to reduce fill-in during factorization
373 @cvar NESTED_DISSECTION: Reorder matrix to improve load balancing during factorization
374 @cvar PASO: PASO solver package
375 @cvar SCSL: SGI SCSL solver library
376 @cvar MKL: Intel's MKL solver library
377 @cvar UMFPACK: the UMFPACK library
378 @cvar TRILINOS: the TRILINOS parallel solver class library from Sandia Natl Labs
379 @cvar ITERATIVE: The default iterative solver
380 @cvar AMG: algebraic multi grid
381 @cvar RILU: recursive ILU
382 @cvar GS: Gauss-Seidel solver
383
384 """
385 DEFAULT= 0
386 DIRECT= 1
387 CHOLEVSKY= 2
388 PCG= 3
389 CR= 4
390 CGS= 5
391 BICGSTAB= 6
392 SSOR= 7
393 ILU0= 8
394 ILUT= 9
395 JACOBI= 10
396 GMRES= 11
397 PRES20= 12
398 LUMPING= 13
399 NO_REORDERING= 17
400 MINIMUM_FILL_IN= 18
401 NESTED_DISSECTION= 19
402 SCSL= 14
403 MKL= 15
404 UMFPACK= 16
405 ITERATIVE= 20
406 PASO= 21
407 AMG= 22
408 RILU = 23
409 TRILINOS = 24
410 NONLINEAR_GMRES = 25
411 TFQMR = 26
412 MINRES = 27
413 GS=28
414
415 SMALL_TOLERANCE=1.e-13
416 PACKAGE_KEY="package"
417 METHOD_KEY="method"
418 SYMMETRY_KEY="symmetric"
419 TOLERANCE_KEY="tolerance"
420 PRECONDITIONER_KEY="preconditioner"
421
422
423 def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):
424 """
425 initializes a linear problem
426
427 @param domain: domain of the PDE
428 @type domain: L{Domain<escript.Domain>}
429 @param numEquations: number of equations. If numEquations==None the number of equations
430 is extracted from the coefficients.
431 @param numSolutions: number of solution components. If numSolutions==None the number of solution components
432 is extracted from the coefficients.
433 @param debug: if True debug informations are printed.
434
435 """
436 super(LinearProblem, self).__init__()
437
438 self.__debug=debug
439 self.__domain=domain
440 self.__numEquations=numEquations
441 self.__numSolutions=numSolutions
442 self.__altered_coefficients=False
443 self.__reduce_equation_order=False
444 self.__reduce_solution_order=False
445 self.__tolerance=1.e-8
446 self.__solver_method=self.DEFAULT
447 self.__solver_package=self.DEFAULT
448 self.__preconditioner=self.DEFAULT
449 self.__is_RHS_valid=False
450 self.__is_operator_valid=False
451 self.__sym=False
452 self.__COEFFICIENTS={}
453 # initialize things:
454 self.resetAllCoefficients()
455 self.__system_type=None
456 self.updateSystemType()
457 # =============================================================================
458 # general stuff:
459 # =============================================================================
460 def __str__(self):
461 """
462 returns string representation of the PDE
463
464 @return: a simple representation of the PDE
465 @rtype: C{str}
466 """
467 return "<LinearProblem %d>"%id(self)
468 # =============================================================================
469 # debug :
470 # =============================================================================
471 def setDebugOn(self):
472 """
473 switches on debugging
474 """
475 self.__debug=not None
476
477 def setDebugOff(self):
478 """
479 switches off debugging
480 """
481 self.__debug=None
482
483 def trace(self,text):
484 """
485 print the text message if debugging is swiched on.
486 @param text: message
487 @type text: C{string}
488 """
489 if self.__debug: print "%s: %s"%(str(self),text)
490
491 # =============================================================================
492 # some service functions:
493 # =============================================================================
494 def introduceCoefficients(self,**coeff):
495 """
496 introduces a new coefficient into the problem.
497
498 use
499
500 self.introduceCoefficients(A=PDECoef(...),B=PDECoef(...))
501
502 to introduce the coefficients M{A} ans M{B}.
503 """
504 for name, type in coeff.items():
505 if not isinstance(type,PDECoef):
506 raise ValueError,"coefficient %s has no type."%name
507 self.__COEFFICIENTS[name]=type
508 self.__COEFFICIENTS[name].resetValue()
509 self.trace("coefficient %s has been introduced."%name)
510
511 def getDomain(self):
512 """
513 returns the domain of the PDE
514
515 @return: the domain of the PDE
516 @rtype: L{Domain<escript.Domain>}
517 """
518 return self.__domain
519
520 def getDim(self):
521 """
522 returns the spatial dimension of the PDE
523
524 @return: the spatial dimension of the PDE domain
525 @rtype: C{int}
526 """
527 return self.getDomain().getDim()
528
529 def getNumEquations(self):
530 """
531 returns the number of equations
532
533 @return: the number of equations
534 @rtype: C{int}
535 @raise UndefinedPDEError: if the number of equations is not be specified yet.
536 """
537 if self.__numEquations==None:
538 raise UndefinedPDEError,"Number of equations is undefined. Please specify argument numEquations."
539 else:
540 return self.__numEquations
541
542 def getNumSolutions(self):
543 """
544 returns the number of unknowns
545
546 @return: the number of unknowns
547 @rtype: C{int}
548 @raise UndefinedPDEError: if the number of unknowns is not be specified yet.
549 """
550 if self.__numSolutions==None:
551 raise UndefinedPDEError,"Number of solution is undefined. Please specify argument numSolutions."
552 else:
553 return self.__numSolutions
554
555 def reduceEquationOrder(self):
556 """
557 return status for order reduction for equation
558
559 @return: return True is reduced interpolation order is used for the represenation of the equation
560 @rtype: L{bool}
561 """
562 return self.__reduce_equation_order
563
564 def reduceSolutionOrder(self):
565 """
566 return status for order reduction for the solution
567
568 @return: return True is reduced interpolation order is used for the represenation of the solution
569 @rtype: L{bool}
570 """
571 return self.__reduce_solution_order
572
573 def getFunctionSpaceForEquation(self):
574 """
575 returns the L{FunctionSpace<escript.FunctionSpace>} used to discretize the equation
576
577 @return: representation space of equation
578 @rtype: L{FunctionSpace<escript.FunctionSpace>}
579 """
580 if self.reduceEquationOrder():
581 return escript.ReducedSolution(self.getDomain())
582 else:
583 return escript.Solution(self.getDomain())
584
585 def getFunctionSpaceForSolution(self):
586 """
587 returns the L{FunctionSpace<escript.FunctionSpace>} used to represent the solution
588
589 @return: representation space of solution
590 @rtype: L{FunctionSpace<escript.FunctionSpace>}
591 """
592 if self.reduceSolutionOrder():
593 return escript.ReducedSolution(self.getDomain())
594 else:
595 return escript.Solution(self.getDomain())
596
597 # =============================================================================
598 # solver settings:
599 # =============================================================================
600 def setSolverMethod(self,solver=None,preconditioner=None):
601 """
602 sets a new solver
603
604 @param solver: sets a new solver method.
605 @type solver: one of L{DEFAULT}, L{ITERATIVE} L{DIRECT}, L{CHOLEVSKY}, L{PCG}, L{CR}, L{CGS}, L{BICGSTAB}, L{SSOR}, L{GMRES}, L{TFQMR}, L{MINRES}, L{PRES20}, L{LUMPING}, L{AMG}
606 @param preconditioner: sets a new solver method.
607 @type preconditioner: one of L{DEFAULT}, L{JACOBI} L{ILU0}, L{ILUT},L{SSOR}, L{RILU}, L{GS}
608
609 @note: the solver method choosen may not be available by the selected package.
610 """
611 if solver==None: solver=self.__solver_method
612 if preconditioner==None: preconditioner=self.__preconditioner
613 if solver==None: solver=self.DEFAULT
614 if preconditioner==None: preconditioner=self.DEFAULT
615 if not (solver,preconditioner)==self.getSolverMethod():
616 self.__solver_method=solver
617 self.__preconditioner=preconditioner
618 self.updateSystemType()
619 self.trace("New solver is %s"%self.getSolverMethodName())
620
621 def getSolverMethodName(self):
622 """
623 returns the name of the solver currently used
624
625 @return: the name of the solver currently used.
626 @rtype: C{string}
627 """
628
629 m=self.getSolverMethod()
630 p=self.getSolverPackage()
631 method=""
632 if m[0]==self.DEFAULT: method="DEFAULT"
633 elif m[0]==self.DIRECT: method= "DIRECT"
634 elif m[0]==self.ITERATIVE: method= "ITERATIVE"
635 elif m[0]==self.CHOLEVSKY: method= "CHOLEVSKY"
636 elif m[0]==self.PCG: method= "PCG"
637 elif m[0]==self.TFQMR: method= "TFQMR"
638 elif m[0]==self.MINRES: method= "MINRES"
639 elif m[0]==self.CR: method= "CR"
640 elif m[0]==self.CGS: method= "CGS"
641 elif m[0]==self.BICGSTAB: method= "BICGSTAB"
642 elif m[0]==self.SSOR: method= "SSOR"
643 elif m[0]==self.GMRES: method= "GMRES"
644 elif m[0]==self.PRES20: method= "PRES20"
645 elif m[0]==self.LUMPING: method= "LUMPING"
646 elif m[0]==self.AMG: method= "AMG"
647 if m[1]==self.DEFAULT: method+="+DEFAULT"
648 elif m[1]==self.JACOBI: method+= "+JACOBI"
649 elif m[1]==self.ILU0: method+= "+ILU0"
650 elif m[1]==self.ILUT: method+= "+ILUT"
651 elif m[1]==self.SSOR: method+= "+SSOR"
652 elif m[1]==self.AMG: method+= "+AMG"
653 elif m[1]==self.RILU: method+= "+RILU"
654 elif m[1]==self.GS: method+= "+GS"
655 if p==self.DEFAULT: package="DEFAULT"
656 elif p==self.PASO: package= "PASO"
657 elif p==self.MKL: package= "MKL"
658 elif p==self.SCSL: package= "SCSL"
659 elif p==self.UMFPACK: package= "UMFPACK"
660 elif p==self.TRILINOS: package= "TRILINOS"
661 else : method="unknown"
662 return "%s solver of %s package"%(method,package)
663
664 def getSolverMethod(self):
665 """
666 returns the solver method and preconditioner used
667
668 @return: the solver method currently be used.
669 @rtype: C{tuple} of C{int}
670 """
671 return self.__solver_method,self.__preconditioner
672
673 def setSolverPackage(self,package=None):
674 """
675 sets a new solver package
676
677 @param package: sets a new solver method.
678 @type package: one of L{DEFAULT}, L{PASO} L{SCSL}, L{MKL}, L{UMFPACK}, L{TRILINOS}
679 """
680 if package==None: package=self.DEFAULT
681 if not package==self.getSolverPackage():
682 self.__solver_package=package
683 self.updateSystemType()
684 self.trace("New solver is %s"%self.getSolverMethodName())
685
686 def getSolverPackage(self):
687 """
688 returns the package of the solver
689
690 @return: the solver package currently being used.
691 @rtype: C{int}
692 """
693 return self.__solver_package
694
695 def isUsingLumping(self):
696 """
697 checks if matrix lumping is used as a solver method
698
699 @return: True is lumping is currently used a solver method.
700 @rtype: C{bool}
701 """
702 return self.getSolverMethod()[0]==self.LUMPING
703
704 def setTolerance(self,tol=1.e-8):
705 """
706 resets the tolerance for the solver method to tol.
707
708 @param tol: new tolerance for the solver. If the tol is lower then the current tolerence
709 the system will be resolved.
710 @type tol: positive C{float}
711 @raise ValueError: if tolerance is not positive.
712 """
713 if not tol>0:
714 raise ValueError,"Tolerance as to be positive"
715 if tol<self.getTolerance(): self.invalidateSolution()
716 self.trace("New tolerance %e"%tol)
717 self.__tolerance=tol
718 return
719
720 def getTolerance(self):
721 """
722 returns the tolerance set for the solution
723
724 @return: tolerance currently used.
725 @rtype: C{float}
726 """
727 return self.__tolerance
728
729 # =============================================================================
730 # symmetry flag:
731 # =============================================================================
732 def isSymmetric(self):
733 """
734 checks if symmetry is indicated.
735
736 @return: True is a symmetric PDE is indicated, otherwise False is returned
737 @rtype: C{bool}
738 """
739 return self.__sym
740
741 def setSymmetryOn(self):
742 """
743 sets the symmetry flag.
744 """
745 if not self.isSymmetric():
746 self.trace("PDE is set to be symmetric")
747 self.__sym=True
748 self.updateSystemType()
749
750 def setSymmetryOff(self):
751 """
752 removes the symmetry flag.
753 """
754 if self.isSymmetric():
755 self.trace("PDE is set to be unsymmetric")
756 self.__sym=False
757 self.updateSystemType()
758
759 def setSymmetryTo(self,flag=False):
760 """
761 sets the symmetry flag to flag
762
763 @param flag: If flag, the symmetry flag is set otherwise the symmetry flag is released.
764 @type flag: C{bool}
765 """
766 if flag:
767 self.setSymmetryOn()
768 else:
769 self.setSymmetryOff()
770
771 # =============================================================================
772 # function space handling for the equation as well as the solution
773 # =============================================================================
774 def setReducedOrderOn(self):
775 """
776 switches on reduced order for solution and equation representation
777
778 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
779 """
780 self.setReducedOrderForSolutionOn()
781 self.setReducedOrderForEquationOn()
782
783 def setReducedOrderOff(self):
784 """
785 switches off reduced order for solution and equation representation
786
787 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
788 """
789 self.setReducedOrderForSolutionOff()
790 self.setReducedOrderForEquationOff()
791
792 def setReducedOrderTo(self,flag=False):
793 """
794 sets order reduction for both solution and equation representation according to flag.
795 @param flag: if flag is True, the order reduction is switched on for both solution and equation representation, otherwise or
796 if flag is not present order reduction is switched off
797 @type flag: C{bool}
798 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
799 """
800 self.setReducedOrderForSolutionTo(flag)
801 self.setReducedOrderForEquationTo(flag)
802
803
804 def setReducedOrderForSolutionOn(self):
805 """
806 switches on reduced order for solution representation
807
808 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
809 """
810 if not self.__reduce_solution_order:
811 if self.__altered_coefficients:
812 raise RuntimeError,"order cannot be altered after coefficients have been defined."
813 self.trace("Reduced order is used to solution representation.")
814 self.__reduce_solution_order=True
815 self.initializeSystem()
816
817 def setReducedOrderForSolutionOff(self):
818 """
819 switches off reduced order for solution representation
820
821 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
822 """
823 if self.__reduce_solution_order:
824 if self.__altered_coefficients:
825 raise RuntimeError,"order cannot be altered after coefficients have been defined."
826 self.trace("Full order is used to interpolate solution.")
827 self.__reduce_solution_order=False
828 self.initializeSystem()
829
830 def setReducedOrderForSolutionTo(self,flag=False):
831 """
832 sets order for test functions according to flag
833
834 @param flag: if flag is True, the order reduction is switched on for solution representation, otherwise or
835 if flag is not present order reduction is switched off
836 @type flag: C{bool}
837 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
838 """
839 if flag:
840 self.setReducedOrderForSolutionOn()
841 else:
842 self.setReducedOrderForSolutionOff()
843
844 def setReducedOrderForEquationOn(self):
845 """
846 switches on reduced order for equation representation
847
848 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
849 """
850 if not self.__reduce_equation_order:
851 if self.__altered_coefficients:
852 raise RuntimeError,"order cannot be altered after coefficients have been defined."
853 self.trace("Reduced order is used for test functions.")
854 self.__reduce_equation_order=True
855 self.initializeSystem()
856
857 def setReducedOrderForEquationOff(self):
858 """
859 switches off reduced order for equation representation
860
861 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
862 """
863 if self.__reduce_equation_order:
864 if self.__altered_coefficients:
865 raise RuntimeError,"order cannot be altered after coefficients have been defined."
866 self.trace("Full order is used for test functions.")
867 self.__reduce_equation_order=False
868 self.initializeSystem()
869
870 def setReducedOrderForEquationTo(self,flag=False):
871 """
872 sets order for test functions according to flag
873
874 @param flag: if flag is True, the order reduction is switched on for equation representation, otherwise or
875 if flag is not present order reduction is switched off
876 @type flag: C{bool}
877 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
878 """
879 if flag:
880 self.setReducedOrderForEquationOn()
881 else:
882 self.setReducedOrderForEquationOff()
883
884 def updateSystemType(self):
885 """
886 reasses the system type and, if a new matrix is needed, resets the system.
887 """
888 new_system_type=self.getRequiredSystemType()
889 if not new_system_type==self.__system_type:
890 self.trace("Matrix type is now %d."%new_system_type)
891 self.__system_type=new_system_type
892 self.initializeSystem()
893 def getSystemType(self):
894 """
895 return the current system type
896 """
897 return self.__system_type
898
899 def checkSymmetricTensor(self,name,verbose=True):
900 """
901 tests a coefficient for symmetry
902
903 @param name: name of the coefficient
904 @type name: C{str}
905 @param verbose: if equal to True or not present a report on coefficients which are breaking the symmetry is printed.
906 @type verbose: C{bool}
907 @return: True if coefficient C{name} is symmetric.
908 @rtype: C{bool}
909 """
910 A=self.getCoefficient(name)
911 verbose=verbose or self.__debug
912 out=True
913 if not A.isEmpty():
914 tol=util.Lsup(A)*self.SMALL_TOLERANCE
915 s=A.getShape()
916 if A.getRank() == 4:
917 if s[0]==s[2] and s[1] == s[3]:
918 for i in range(s[0]):
919 for j in range(s[1]):
920 for k in range(s[2]):
921 for l in range(s[3]):
922 if util.Lsup(A[i,j,k,l]-A[k,l,i,j])>tol:
923 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)
924 out=False
925 else:
926 if verbose: print "non-symmetric problem because of inappropriate shape %s of coefficient %s."%(s,name)
927 out=False
928 elif A.getRank() == 2:
929 if s[0]==s[1]:
930 for j in range(s[0]):
931 for l in range(s[1]):
932 if util.Lsup(A[j,l]-A[l,j])>tol:
933 if verbose: print "non-symmetric problem because %s[%d,%d]!=%s[%d,%d]"%(name,j,l,name,l,j)
934 out=False
935 else:
936 if verbose: print "non-symmetric problem because of inappropriate shape %s of coefficient %s."%(s,name)
937 out=False
938 elif A.getRank() == 0:
939 pass
940 else:
941 raise ValueError,"Cannot check rank %s of %s."%(A.getRank(),name)
942 return out
943
944 def checkReciprocalSymmetry(self,name0,name1,verbose=True):
945 """
946 test two coefficients for resciprocal symmetry
947
948 @param name0: name of the first coefficient
949 @type name: C{str}
950 @param name1: name of the second coefficient
951 @type name: C{str}
952 @param verbose: if equal to True or not present a report on coefficients which are breaking the symmetry is printed.
953 @type verbose: C{bool}
954 @return: True if coefficients C{name0} and C{name1} are resciprocally symmetric.
955 @rtype: C{bool}
956 """
957 B=self.getCoefficient(name0)
958 C=self.getCoefficient(name1)
959 verbose=verbose or self.__debug
960 out=True
961 if B.isEmpty() and not C.isEmpty():
962 if verbose: print "non-symmetric problem because %s is not present but %s is"%(name0,name1)
963 out=False
964 elif not B.isEmpty() and C.isEmpty():
965 if verbose: print "non-symmetric problem because %s is not present but %s is"%(name0,name1)
966 out=False
967 elif not B.isEmpty() and not C.isEmpty():
968 sB=B.getShape()
969 sC=C.getShape()
970 tol=(util.Lsup(B)+util.Lsup(C))*self.SMALL_TOLERANCE/2.
971 if len(sB) != len(sC):
972 if verbose: print "non-symmetric problem because ranks of %s (=%s) and %s (=%s) are different."%(name0,len(sB),name1,len(sC))
973 out=False
974 else:
975 if len(sB)==0:
976 if util.Lsup(B-C)>tol:
977 if verbose: print "non-symmetric problem because %s!=%s"%(name0,name1)
978 out=False
979 elif len(sB)==1:
980 if sB[0]==sC[0]:
981 for j in range(sB[0]):
982 if util.Lsup(B[j]-C[j])>tol:
983 if verbose: print "non-symmetric PDE because %s[%d]!=%s[%d]"%(name0,j,name1,j)
984 out=False
985 else:
986 if verbose: print "non-symmetric problem because of inappropriate shapes %s and %s of coefficients %s and %s, respectively."%(sB,sC,name0,name1)
987 elif len(sB)==3:
988 if sB[0]==sC[1] and sB[1]==sC[2] and sB[2]==sC[0]:
989 for i in range(sB[0]):
990 for j in range(sB[1]):
991 for k in range(sB[2]):
992 if util.Lsup(B[i,j,k]-C[k,i,j])>tol:
993 if verbose: print "non-symmetric problem because %s[%d,%d,%d]!=%s[%d,%d,%d]"%(name2,i,j,k,name1,k,i,j)
994 out=False
995 else:
996 if verbose: print "non-symmetric problem because of inappropriate shapes %s and %s of coefficients %s and %s, respectively."%(sB,sC,name0,name1)
997 else:
998 raise ValueError,"Cannot check rank %s of %s and %s."%(len(sB),name0,name1)
999 return out
1000
1001 def getCoefficient(self,name):
1002 """
1003 returns the value of the coefficient name
1004
1005 @param name: name of the coefficient requested.
1006 @type name: C{string}
1007 @return: the value of the coefficient name
1008 @rtype: L{Data<escript.Data>}
1009 @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1010 """
1011 if self.hasCoefficient(name):
1012 return self.__COEFFICIENTS[name].getValue()
1013 else:
1014 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1015
1016 def hasCoefficient(self,name):
1017 """
1018 return True if name is the name of a coefficient
1019
1020 @param name: name of the coefficient enquired.
1021 @type name: C{string}
1022 @return: True if name is the name of a coefficient of the general PDE. Otherwise False.
1023 @rtype: C{bool}
1024 """
1025 return self.__COEFFICIENTS.has_key(name)
1026
1027 def createCoefficient(self, name):
1028 """
1029 create a L{Data<escript.Data>} object corresponding to coefficient name
1030
1031 @return: a coefficient name initialized to 0.
1032 @rtype: L{Data<escript.Data>}
1033 @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1034 """
1035 if self.hasCoefficient(name):
1036 return escript.Data(0.,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name))
1037 else:
1038 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1039
1040 def getFunctionSpaceForCoefficient(self,name):
1041 """
1042 return the L{FunctionSpace<escript.FunctionSpace>} to be used for coefficient name
1043
1044 @param name: name of the coefficient enquired.
1045 @type name: C{string}
1046 @return: the function space to be used for coefficient name
1047 @rtype: L{FunctionSpace<escript.FunctionSpace>}
1048 @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1049 """
1050 if self.hasCoefficient(name):
1051 return self.__COEFFICIENTS[name].getFunctionSpace(self.getDomain())
1052 else:
1053 raise ValueError,"unknown coefficient %s requested"%name
1054
1055 def getShapeOfCoefficient(self,name):
1056 """
1057 return the shape of the coefficient name
1058
1059 @param name: name of the coefficient enquired.
1060 @type name: C{string}
1061 @return: the shape of the coefficient name
1062 @rtype: C{tuple} of C{int}
1063 @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1064 """
1065 if self.hasCoefficient(name):
1066 return self.__COEFFICIENTS[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())
1067 else:
1068 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1069
1070 def resetAllCoefficients(self):
1071 """
1072 resets all coefficients to there default values.
1073 """
1074 for i in self.__COEFFICIENTS.iterkeys():
1075 self.__COEFFICIENTS[i].resetValue()
1076
1077 def alteredCoefficient(self,name):
1078 """
1079 announce that coefficient name has been changed
1080
1081 @param name: name of the coefficient enquired.
1082 @type name: C{string}
1083 @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1084 @note: if name is q or r, the method will not trigger a rebuilt of the system as constraints are applied to the solved system.
1085 """
1086 if self.hasCoefficient(name):
1087 self.trace("Coefficient %s has been altered."%name)
1088 if not ((name=="q" or name=="r") and self.isUsingLumping()):
1089 if self.__COEFFICIENTS[name].isAlteringOperator(): self.invalidateOperator()
1090 if self.__COEFFICIENTS[name].isAlteringRightHandSide(): self.invalidateRightHandSide()
1091 else:
1092 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1093
1094 def validSolution(self):
1095 """
1096 markes the solution as valid
1097 """
1098 self.__is_solution_valid=True
1099
1100 def invalidateSolution(self):
1101 """
1102 indicates the PDE has to be resolved if the solution is requested
1103 """
1104 self.trace("System will be resolved.")
1105 self.__is_solution_valid=False
1106
1107 def isSolutionValid(self):
1108 """
1109 returns True is the solution is still valid.
1110 """
1111 return self.__is_solution_valid
1112
1113 def validOperator(self):
1114 """
1115 markes the operator as valid
1116 """
1117 self.__is_operator_valid=True
1118
1119 def invalidateOperator(self):
1120 """
1121 indicates the operator has to be rebuilt next time it is used
1122 """
1123 self.trace("Operator will be rebuilt.")
1124 self.invalidateSolution()
1125 self.__is_operator_valid=False
1126
1127 def isOperatorValid(self):
1128 """
1129 returns True is the operator is still valid.
1130 """
1131 return self.__is_operator_valid
1132
1133 def validRightHandSide(self):
1134 """
1135 markes the right hand side as valid
1136 """
1137 self.__is_RHS_valid=True
1138
1139 def invalidateRightHandSide(self):
1140 """
1141 indicates the right hand side has to be rebuild next time it is used
1142 """
1143 if self.isRightHandSideValid(): self.trace("Right hand side has to be rebuilt.")
1144 self.invalidateSolution()
1145 self.__is_RHS_valid=False
1146
1147 def isRightHandSideValid(self):
1148 """
1149 returns True is the operator is still valid.
1150 """
1151 return self.__is_RHS_valid
1152
1153 def invalidateSystem(self):
1154 """
1155 annonced that everthing has to be rebuild:
1156 """
1157 if self.isRightHandSideValid(): self.trace("System has to be rebuilt.")
1158 self.invalidateSolution()
1159 self.invalidateOperator()
1160 self.invalidateRightHandSide()
1161
1162 def isSystemValid(self):
1163 """
1164 returns true if the system (including solution) is still vaild:
1165 """
1166 return self.isSolutionValid() and self.isOperatorValid() and self.isRightHandSideValid()
1167
1168 def initializeSystem(self):
1169 """
1170 resets the system
1171 """
1172 self.trace("New System has been created.")
1173 self.__operator=escript.Operator()
1174 self.__righthandside=escript.Data()
1175 self.__solution=escript.Data()
1176 self.invalidateSystem()
1177
1178 def getOperator(self):
1179 """
1180 returns the operator of the linear problem
1181
1182 @return: the operator of the problem
1183 """
1184 return self.getSystem()[0]
1185
1186 def getRightHandSide(self):
1187 """
1188 returns the right hand side of the linear problem
1189
1190 @return: the right hand side of the problem
1191 @rtype: L{Data<escript.Data>}
1192 """
1193 return self.getSystem()[1]
1194
1195 def createRightHandSide(self):
1196 """
1197 returns an instance of a new right hand side
1198 """
1199 self.trace("New right hand side is allocated.")
1200 if self.getNumEquations()>1:
1201 return escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True)
1202 else:
1203 return escript.Data(0.,(),self.getFunctionSpaceForEquation(),True)
1204
1205 def createSolution(self):
1206 """
1207 returns an instance of a new solution
1208 """
1209 self.trace("New solution is allocated.")
1210 if self.getNumSolutions()>1:
1211 return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)
1212 else:
1213 return escript.Data(0.,(),self.getFunctionSpaceForSolution(),True)
1214
1215 def resetSolution(self):
1216 """
1217 sets the solution to zero
1218 """
1219 if self.__solution.isEmpty():
1220 self.__solution=self.createSolution()
1221 else:
1222 self.__solution.setToZero()
1223 self.trace("Solution is reset to zero.")
1224 def setSolution(self,u):
1225 """
1226 sets the solution assuming that makes the sytem valid.
1227 """
1228 self.__solution=u
1229 self.validSolution()
1230
1231 def getCurrentSolution(self):
1232 """
1233 returns the solution in its current state
1234 """
1235 if self.__solution.isEmpty(): self.__solution=self.createSolution()
1236 return self.__solution
1237
1238 def resetRightHandSide(self):
1239 """
1240 sets the right hand side to zero
1241 """
1242 if self.__righthandside.isEmpty():
1243 self.__righthandside=self.createRightHandSide()
1244 else:
1245 self.__righthandside.setToZero()
1246 self.trace("Right hand side is reset to zero.")
1247
1248 def getCurrentRightHandSide(self):
1249 """
1250 returns the right hand side in its current state
1251 """
1252 if self.__righthandside.isEmpty(): self.__righthandside=self.createRightHandSide()
1253 return self.__righthandside
1254
1255
1256 def resetOperator(self):
1257 """
1258 makes sure that the operator is instantiated and returns it initialized by zeros
1259 """
1260 if self.__operator.isEmpty():
1261 if self.isUsingLumping():
1262 self.__operator=self.createSolution()
1263 else:
1264 self.__operator=self.createOperator()
1265 else:
1266 if self.isUsingLumping():
1267 self.__operator.setToZero()
1268 else:
1269 self.__operator.resetValues()
1270 self.trace("Operator reset to zero")
1271
1272 def getCurrentOperator(self):
1273 """
1274 returns the operator in its current state
1275 """
1276 if self.__operator.isEmpty(): self.resetOperator()
1277 return self.__operator
1278
1279 def setValue(self,**coefficients):
1280 """
1281 sets new values to coefficients
1282
1283 @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1284 """
1285 # check if the coefficients are legal:
1286 for i in coefficients.iterkeys():
1287 if not self.hasCoefficient(i):
1288 raise IllegalCoefficient,"Attempt to set unknown coefficient %s"%i
1289 # if the number of unknowns or equations is still unknown we try to estimate them:
1290 if self.__numEquations==None or self.__numSolutions==None:
1291 for i,d in coefficients.iteritems():
1292 if hasattr(d,"shape"):
1293 s=d.shape
1294 elif hasattr(d,"getShape"):
1295 s=d.getShape()
1296 else:
1297 s=numarray.array(d).shape
1298 if s!=None:
1299 # get number of equations and number of unknowns:
1300 res=self.__COEFFICIENTS[i].estimateNumEquationsAndNumSolutions(self.getDomain(),s)
1301 if res==None:
1302 raise IllegalCoefficientValue,"Illegal shape %s of coefficient %s"%(s,i)
1303 else:
1304 if self.__numEquations==None: self.__numEquations=res[0]
1305 if self.__numSolutions==None: self.__numSolutions=res[1]
1306 if self.__numEquations==None: raise UndefinedPDEError,"unidententified number of equations"
1307 if self.__numSolutions==None: raise UndefinedPDEError,"unidententified number of solutions"
1308 # now we check the shape of the coefficient if numEquations and numSolutions are set:
1309 for i,d in coefficients.iteritems():
1310 try:
1311 self.__COEFFICIENTS[i].setValue(self.getDomain(),
1312 self.getNumEquations(),self.getNumSolutions(),
1313 self.reduceEquationOrder(),self.reduceSolutionOrder(),d)
1314 self.alteredCoefficient(i)
1315 except IllegalCoefficientFunctionSpace,m:
1316 # if the function space is wrong then we try the reduced version:
1317 i_red=i+"_reduced"
1318 if (not i_red in coefficients.keys()) and i_red in self.__COEFFICIENTS.keys():
1319 try:
1320 self.__COEFFICIENTS[i_red].setValue(self.getDomain(),
1321 self.getNumEquations(),self.getNumSolutions(),
1322 self.reduceEquationOrder(),self.reduceSolutionOrder(),d)
1323 self.alteredCoefficient(i_red)
1324 except IllegalCoefficientValue,m:
1325 raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
1326 except IllegalCoefficientFunctionSpace,m:
1327 raise IllegalCoefficientFunctionSpace("Coefficient %s:%s"%(i,m))
1328 else:
1329 raise IllegalCoefficientFunctionSpace("Coefficient %s:%s"%(i,m))
1330 except IllegalCoefficientValue,m:
1331 raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
1332 self.__altered_coefficients=True
1333
1334 # ====================================================================================
1335 # methods that are typically overwritten when implementing a particular linear problem
1336 # ====================================================================================
1337 def getRequiredSystemType(self):
1338 """
1339 returns the system type which needs to be used by the current set up.
1340
1341 @note: Typically this method is overwritten when implementing a particular linear problem.
1342 """
1343 return 0
1344
1345 def createOperator(self):
1346 """
1347 returns an instance of a new operator
1348
1349 @note: This method is overwritten when implementing a particular linear problem.
1350 """
1351 return escript.Operator()
1352
1353 def checkSymmetry(self,verbose=True):
1354 """
1355 test the PDE for symmetry.
1356
1357 @param verbose: if equal to True or not present a report on coefficients which are breaking the symmetry is printed.
1358 @type verbose: C{bool}
1359 @return: True if the problem is symmetric.
1360 @rtype: C{bool}
1361 @note: Typically this method is overwritten when implementing a particular linear problem.
1362 """
1363 out=True
1364 return out
1365
1366 def getSolution(self,**options):
1367 """
1368 returns the solution of the problem
1369 @return: the solution
1370 @rtype: L{Data<escript.Data>}
1371
1372 @note: This method is overwritten when implementing a particular linear problem.
1373 """
1374 return self.getCurrentSolution()
1375
1376 def getSystem(self):
1377 """
1378 return the operator and right hand side of the PDE
1379
1380 @return: the discrete version of the PDE
1381 @rtype: C{tuple} of L{Operator,<escript.Operator>} and L{Data<escript.Data>}.
1382
1383 @note: This method is overwritten when implementing a particular linear problem.
1384 """
1385 return (self.getCurrentOperator(), self.getCurrentRightHandSide())
1386 #=====================
1387
1388 class LinearPDE(LinearProblem):
1389 """
1390 This class is used to define a general linear, steady, second order PDE
1391 for an unknown function M{u} on a given domain defined through a L{Domain<escript.Domain>} object.
1392
1393 For a single PDE with a solution with a single component the linear PDE is defined in the following form:
1394
1395 M{-(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)}
1396
1397
1398 where M{grad(F)} denotes the spatial derivative of M{F}. Einstein's summation convention,
1399 ie. summation over indexes appearing twice in a term of a sum is performed, is used.
1400 The coefficients M{A}, M{B}, M{C}, M{D}, M{X} and M{Y} have to be specified through L{Data<escript.Data>} objects in the
1401 L{Function<escript.Function>} and the coefficients M{A_reduced}, M{B_reduced}, M{C_reduced}, M{D_reduced}, M{X_reduced} and M{Y_reduced} have to be specified through L{Data<escript.Data>} objects in the
1402 L{ReducedFunction<escript.ReducedFunction>}. It is also allowd to use objects that can be converted into
1403 such L{Data<escript.Data>} objects. M{A} and M{A_reduced} are rank two, M{B}, M{C}, M{X},
1404 M{B_reduced}, M{C_reduced} and M{X_reduced} are rank one and M{D}, M{D_reduced} M{Y} and M{Y_reduced} are scalar.
1405
1406 The following natural boundary conditions are considered:
1407
1408 M{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}
1409
1410 where M{n} is the outer normal field. Notice that the coefficients M{A}, M{A_reduced}, M{B}, M{B_reduced}, M{X} and M{X_reduced} are defined in the PDE. The coefficients M{d} and M{y} and are each a scalar in the L{FunctionOnBoundary<escript.FunctionOnBoundary>} and the coefficients M{d_reduced} and M{y_reduced} and are each a scalar in the L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
1411
1412
1413 Constraints for the solution prescribing the value of the solution at certain locations in the domain. They have the form
1414
1415 M{u=r} where M{q>0}
1416
1417 M{r} and M{q} are each scalar where M{q} is the characteristic function defining where the constraint is applied.
1418 The constraints override any other condition set by the PDE or the boundary condition.
1419
1420 The PDE is symmetrical if
1421
1422 M{A[i,j]=A[j,i]} and M{B[j]=C[j]} and M{A_reduced[i,j]=A_reduced[j,i]} and M{B_reduced[j]=C_reduced[j]}
1423
1424 For a system of PDEs and a solution with several components the PDE has the form
1425
1426 M{-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] }
1427
1428 M{A} and M{A_reduced} are of rank four, M{B}, M{B_reduced}, M{C} and M{C_reduced} are each of rank three, M{D}, M{D_reduced}, M{X_reduced} and M{X} are each a rank two and M{Y} and M{Y_reduced} are of rank one.
1429 The natural boundary conditions take the form:
1430
1431 M{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]}
1432
1433
1434 The coefficient M{d} is a rank two and M{y} is a rank one both in the L{FunctionOnBoundary<escript.FunctionOnBoundary>}. Constraints take the form and the coefficients M{d_reduced} is a rank two and M{y_reduced} is a rank one both in the L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
1435
1436 Constraints take the form
1437
1438 M{u[i]=r[i]} where M{q[i]>0}
1439
1440 M{r} and M{q} are each rank one. Notice that at some locations not necessarily all components must have a constraint.
1441
1442 The system of PDEs is symmetrical if
1443
1444 - M{A[i,j,k,l]=A[k,l,i,j]}
1445 - M{A_reduced[i,j,k,l]=A_reduced[k,l,i,j]}
1446 - M{B[i,j,k]=C[k,i,j]}
1447 - M{B_reduced[i,j,k]=C_reduced[k,i,j]}
1448 - M{D[i,k]=D[i,k]}
1449 - M{D_reduced[i,k]=D_reduced[i,k]}
1450 - M{d[i,k]=d[k,i]}
1451 - M{d_reduced[i,k]=d_reduced[k,i]}
1452
1453 L{LinearPDE} also supports solution discontinuities over a contact region in the domain. To specify the conditions across the
1454 discontinuity we are using the generalised flux M{J} which is in the case of a systems of PDEs and several components of the solution
1455 defined as
1456
1457 M{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]}
1458
1459 For the case of single solution component and single PDE M{J} is defined
1460
1461 M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[j]+(B[i]+B_reduced[i])*u-X[i]-X_reduced[i]}
1462
1463 In the context of discontinuities M{n} denotes the normal on the discontinuity pointing from side 0 towards side 1
1464 calculated from L{getNormal<escript.FunctionSpace.getNormal>} of L{FunctionOnContactZero<escript.FunctionOnContactZero>}. For a system of PDEs
1465 the contact condition takes the form
1466
1467 M{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]}
1468
1469 where M{J0} and M{J1} are the fluxes on side 0 and side 1 of the discontinuity, respectively. M{jump(u)}, which is the difference
1470 of the solution at side 1 and at side 0, denotes the jump of M{u} across discontinuity along the normal calculated by
1471 L{jump<util.jump>}.
1472 The coefficient M{d_contact} is a rank two and M{y_contact} is a rank one both in the L{FunctionOnContactZero<escript.FunctionOnContactZero>} or L{FunctionOnContactOne<escript.FunctionOnContactOne>}.
1473 The coefficient M{d_contact_reduced} is a rank two and M{y_contact_reduced} is a rank one both in the L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}.
1474 In case of a single PDE and a single component solution the contact condition takes the form
1475
1476 M{n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)}
1477
1478 In this case the coefficient M{d_contact} and M{y_contact} are each scalar both in the L{FunctionOnContactZero<escript.FunctionOnContactZero>} or L{FunctionOnContactOne<escript.FunctionOnContactOne>} and the coefficient M{d_contact_reduced} and M{y_contact_reduced} are each scalar both in the L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}
1479
1480 typical usage:
1481
1482 p=LinearPDE(dom)
1483 p.setValue(A=kronecker(dom),D=1,Y=0.5)
1484 u=p.getSolution()
1485
1486 """
1487
1488
1489 def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):
1490 """
1491 initializes a new linear PDE
1492
1493 @param domain: domain of the PDE
1494 @type domain: L{Domain<escript.Domain>}
1495 @param numEquations: number of equations. If numEquations==None the number of equations
1496 is extracted from the PDE coefficients.
1497 @param numSolutions: number of solution components. If numSolutions==None the number of solution components
1498 is extracted from the PDE coefficients.
1499 @param debug: if True debug informations are printed.
1500
1501 """
1502 super(LinearPDE, self).__init__(domain,numEquations,numSolutions,debug)
1503 #
1504 # the coefficients of the PDE:
1505 #
1506 self.introduceCoefficients(
1507 A=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
1508 B=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1509 C=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
1510 D=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1511 X=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
1512 Y=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1513 d=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1514 y=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1515 d_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1516 y_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1517 A_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
1518 B_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1519 C_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
1520 D_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1521 X_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
1522 Y_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1523 d_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1524 y_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1525 d_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1526 y_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1527 r=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.RIGHTHANDSIDE),
1528 q=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.BOTH) )
1529
1530 def __str__(self):
1531 """
1532 returns string representation of the PDE
1533
1534 @return: a simple representation of the PDE
1535 @rtype: C{str}
1536 """
1537 return "<LinearPDE %d>"%id(self)
1538
1539 def getRequiredSystemType(self):
1540 """
1541 returns the system type which needs to be used by the current set up.
1542 """
1543 return self.getDomain().getSystemMatrixTypeId(self.getSolverMethod()[0],self.getSolverPackage(),self.isSymmetric())
1544
1545 def checkSymmetry(self,verbose=True):
1546 """
1547 test the PDE for symmetry.
1548
1549 @param verbose: if equal to True or not present a report on coefficients which are breaking the symmetry is printed.
1550 @type verbose: C{bool}
1551 @return: True if the PDE is symmetric.
1552 @rtype: L{bool}
1553 @note: This is a very expensive operation. It should be used for degugging only! The symmetry flag is not altered.
1554 """
1555 out=True
1556 out=out and self.checkSymmetricTensor("A", verbose)
1557 out=out and self.checkSymmetricTensor("A_reduced", verbose)
1558 out=out and self.checkReciprocalSymmetry("B","C", verbose)
1559 out=out and self.checkReciprocalSymmetry("B_reduced","C_reduced", verbose)
1560 out=out and self.checkSymmetricTensor("D", verbose)
1561 out=out and self.checkSymmetricTensor("D_reduced", verbose)
1562 out=out and self.checkSymmetricTensor("d", verbose)
1563 out=out and self.checkSymmetricTensor("d_reduced", verbose)
1564 out=out and self.checkSymmetricTensor("d_contact", verbose)
1565 out=out and self.checkSymmetricTensor("d_contact_reduced", verbose)
1566 return out
1567
1568 def createOperator(self):
1569 """
1570 returns an instance of a new operator
1571 """
1572 self.trace("New operator is allocated.")
1573 return self.getDomain().newOperator( \
1574 self.getNumEquations(), \
1575 self.getFunctionSpaceForEquation(), \
1576 self.getNumSolutions(), \
1577 self.getFunctionSpaceForSolution(), \
1578 self.getSystemType())
1579
1580 def getSolution(self,**options):
1581 """
1582 returns the solution of the PDE
1583
1584 @return: the solution
1585 @rtype: L{Data<escript.Data>}
1586 @param options: solver options
1587 @keyword verbose: True to get some information during PDE solution
1588 @type verbose: C{bool}
1589 @keyword reordering: reordering scheme to be used during elimination. Allowed values are
1590 L{NO_REORDERING}, L{MINIMUM_FILL_IN}, L{NESTED_DISSECTION}
1591 @keyword iter_max: maximum number of iteration steps allowed.
1592 @keyword drop_tolerance: threshold for drupping in L{ILUT}
1593 @keyword drop_storage: maximum of allowed memory in L{ILUT}
1594 @keyword truncation: maximum number of residuals in L{GMRES}
1595 @keyword restart: restart cycle length in L{GMRES}
1596 """
1597 if not self.isSolutionValid():
1598 mat,f=self.getSystem()
1599 if self.isUsingLumping():
1600 self.setSolution(f*1/mat)
1601 else:
1602 options[self.TOLERANCE_KEY]=self.getTolerance()
1603 options[self.METHOD_KEY]=self.getSolverMethod()[0]
1604 options[self.PRECONDITIONER_KEY]=self.getSolverMethod()[1]
1605 options[self.PACKAGE_KEY]=self.getSolverPackage()
1606 options[self.SYMMETRY_KEY]=self.isSymmetric()
1607 self.trace("PDE is resolved.")
1608 self.trace("solver options: %s"%str(options))
1609 self.setSolution(mat.solve(f,options))
1610 return self.getCurrentSolution()
1611
1612 def getSystem(self):
1613 """
1614 return the operator and right hand side of the PDE
1615
1616 @return: the discrete version of the PDE
1617 @rtype: C{tuple} of L{Operator,<escript.Operator>} and L{Data<escript.Data>}.
1618 """
1619 if not self.isOperatorValid() or not self.isRightHandSideValid():
1620 if self.isUsingLumping():
1621 if not self.isOperatorValid():
1622 if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():
1623 raise TypeError,"Lumped matrix requires same order for equations and unknowns"
1624 if not self.getCoefficient("A").isEmpty():
1625 raise ValueError,"coefficient A in lumped matrix may not be present."
1626 if not self.getCoefficient("B").isEmpty():
1627 raise ValueError,"coefficient B in lumped matrix may not be present."
1628 if not self.getCoefficient("C").isEmpty():
1629 raise ValueError,"coefficient C in lumped matrix may not be present."
1630 if not self.getCoefficient("d_contact").isEmpty():
1631 raise ValueError,"coefficient d_contact in lumped matrix may not be present."
1632 if not self.getCoefficient("A_reduced").isEmpty():
1633 raise ValueError,"coefficient A_reduced in lumped matrix may not be present."
1634 if not self.getCoefficient("B_reduced").isEmpty():
1635 raise ValueError,"coefficient B_reduced in lumped matrix may not be present."
1636 if not self.getCoefficient("C_reduced").isEmpty():
1637 raise ValueError,"coefficient C_reduced in lumped matrix may not be present."
1638 if not self.getCoefficient("d_contact_reduced").isEmpty():
1639 raise ValueError,"coefficient d_contact_reduced in lumped matrix may not be present."
1640 D=self.getCoefficient("D")
1641 d=self.getCoefficient("d")
1642 D_reduced=self.getCoefficient("D_reduced")
1643 d_reduced=self.getCoefficient("d_reduced")
1644 if not D.isEmpty():
1645 if self.getNumSolutions()>1:
1646 D_times_e=util.matrix_mult(D,numarray.ones((self.getNumSolutions(),)))
1647 else:
1648 D_times_e=D
1649 else:
1650 D_times_e=escript.Data()
1651 if not d.isEmpty():
1652 if self.getNumSolutions()>1:
1653 d_times_e=util.matrix_mult(d,numarray.ones((self.getNumSolutions(),)))
1654 else:
1655 d_times_e=d
1656 else:
1657 d_times_e=escript.Data()
1658
1659 if not D_reduced.isEmpty():
1660 if self.getNumSolutions()>1:
1661 D_reduced_times_e=util.matrix_mult(D_reduced,numarray.ones((self.getNumSolutions(),)))
1662 else:
1663 D_reduced_times_e=D_reduced
1664 else:
1665 D_reduced_times_e=escript.Data()
1666 if not d_reduced.isEmpty():
1667 if self.getNumSolutions()>1:
1668 d_reduced_times_e=util.matrix_mult(d_reduced,numarray.ones((self.getNumSolutions(),)))
1669 else:
1670 d_reduced_times_e=d_reduced
1671 else:
1672 d_reduced_times_e=escript.Data()
1673
1674 self.resetOperator()
1675 operator=self.getCurrentOperator()
1676 if False and hasattr(self.getDomain(), "addPDEToLumpedSystem") :
1677 self.getDomain().addPDEToLumpedSystem(operator, D_times_e, d_times_e)
1678 self.getDomain().addPDEToLumpedSystem(operator, D_reduced_times_e, d_reduced_times_e)
1679 else:
1680 self.getDomain().addPDEToRHS(operator, \
1681 escript.Data(), \
1682 D_times_e, \
1683 d_times_e,\
1684 escript.Data())
1685 self.getDomain().addPDEToRHS(operator, \
1686 escript.Data(), \
1687 D_reduced_times_e, \
1688 d_reduced_times_e,\
1689 escript.Data())
1690 self.trace("New lumped operator has been built.")
1691 if not self.isRightHandSideValid():
1692 self.resetRightHandSide()
1693 righthandside=self.getCurrentRightHandSide()
1694 self.getDomain().addPDEToRHS(righthandside, \
1695 self.getCoefficient("X"), \
1696 self.getCoefficient("Y"),\
1697 self.getCoefficient("y"),\
1698 self.getCoefficient("y_contact"))
1699 self.getDomain().addPDEToRHS(righthandside, \
1700 self.getCoefficient("X_reduced"), \
1701 self.getCoefficient("Y_reduced"),\
1702 self.getCoefficient("y_reduced"),\
1703 self.getCoefficient("y_contact_reduced"))
1704 self.trace("New right hand side as been built.")
1705 self.validRightHandSide()
1706 self.insertConstraint(rhs_only=False)
1707 self.validOperator()
1708 else:
1709 if not self.isOperatorValid() and not self.isRightHandSideValid():
1710 self.resetRightHandSide()
1711 righthandside=self.getCurrentRightHandSide()
1712 self.resetOperator()
1713 operator=self.getCurrentOperator()
1714 self.getDomain().addPDEToSystem(operator,righthandside, \
1715 self.getCoefficient("A"), \
1716 self.getCoefficient("B"), \
1717 self.getCoefficient("C"), \
1718 self.getCoefficient("D"), \
1719 self.getCoefficient("X"), \
1720 self.getCoefficient("Y"), \
1721 self.getCoefficient("d"), \
1722 self.getCoefficient("y"), \
1723 self.getCoefficient("d_contact"), \
1724 self.getCoefficient("y_contact"))
1725 self.getDomain().addPDEToSystem(operator,righthandside, \
1726 self.getCoefficient("A_reduced"), \
1727 self.getCoefficient("B_reduced"), \
1728 self.getCoefficient("C_reduced"), \
1729 self.getCoefficient("D_reduced"), \
1730 self.getCoefficient("X_reduced"), \
1731 self.getCoefficient("Y_reduced"), \
1732 self.getCoefficient("d_reduced"), \
1733 self.getCoefficient("y_reduced"), \
1734 self.getCoefficient("d_contact_reduced"), \
1735 self.getCoefficient("y_contact_reduced"))
1736 self.insertConstraint(rhs_only=False)
1737 self.trace("New system has been built.")
1738 self.validOperator()
1739 self.validRightHandSide()
1740 elif not self.isRightHandSideValid():
1741 self.resetRightHandSide()
1742 righthandside=self.getCurrentRightHandSide()
1743 self.getDomain().addPDEToRHS(righthandside,
1744 self.getCoefficient("X"), \
1745 self.getCoefficient("Y"),\
1746 self.getCoefficient("y"),\
1747 self.getCoefficient("y_contact"))
1748 self.getDomain().addPDEToRHS(righthandside,
1749 self.getCoefficient("X_reduced"), \
1750 self.getCoefficient("Y_reduced"),\
1751 self.getCoefficient("y_reduced"),\
1752 self.getCoefficient("y_contact_reduced"))
1753 self.insertConstraint(rhs_only=True)
1754 self.trace("New right hand side has been built.")
1755 self.validRightHandSide()
1756 elif not self.isOperatorValid():
1757 self.resetOperator()
1758 operator=self.getCurrentOperator()
1759 self.getDomain().addPDEToSystem(operator,escript.Data(), \
1760 self.getCoefficient("A"), \
1761 self.getCoefficient("B"), \
1762 self.getCoefficient("C"), \
1763 self.getCoefficient("D"), \
1764 escript.Data(), \
1765 escript.Data(), \
1766 self.getCoefficient("d"), \
1767 escript.Data(),\
1768 self.getCoefficient("d_contact"), \
1769 escript.Data())
1770 self.getDomain().addPDEToSystem(operator,escript.Data(), \
1771 self.getCoefficient("A_reduced"), \
1772 self.getCoefficient("B_reduced"), \
1773 self.getCoefficient("C_reduced"), \
1774 self.getCoefficient("D_reduced"), \
1775 escript.Data(), \
1776 escript.Data(), \
1777 self.getCoefficient("d_reduced"), \
1778 escript.Data(),\
1779 self.getCoefficient("d_contact_reduced"), \
1780 escript.Data())
1781 self.insertConstraint(rhs_only=False)
1782 self.trace("New operator has been built.")
1783 self.validOperator()
1784 return (self.getCurrentOperator(), self.getCurrentRightHandSide())
1785
1786 def insertConstraint(self, rhs_only=False):
1787 """
1788 applies the constraints defined by q and r to the PDE
1789
1790 @param rhs_only: if True only the right hand side is altered by the constraint.
1791 @type rhs_only: C{bool}
1792 """
1793 q=self.getCoefficient("q")
1794 r=self.getCoefficient("r")
1795 righthandside=self.getCurrentRightHandSide()
1796 operator=self.getCurrentOperator()
1797
1798 if not q.isEmpty():
1799 if r.isEmpty():
1800 r_s=self.createSolution()
1801 else:
1802 r_s=r
1803 if not rhs_only and not operator.isEmpty():
1804 if self.isUsingLumping():
1805 operator.copyWithMask(escript.Data(1.,q.getShape(),q.getFunctionSpace()),q)
1806 else:
1807 row_q=escript.Data(q,self.getFunctionSpaceForEquation())
1808 col_q=escript.Data(q,self.getFunctionSpaceForSolution())
1809 u=self.createSolution()
1810 u.copyWithMask(r_s,col_q)
1811 righthandside-=operator*u
1812 operator.nullifyRowsAndCols(row_q,col_q,1.)
1813 righthandside.copyWithMask(r_s,q)
1814
1815 def setValue(self,**coefficients):
1816 """
1817 sets new values to coefficients
1818
1819 @param coefficients: new values assigned to coefficients
1820 @keyword A: value for coefficient A.
1821 @type A: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1822 @keyword A_reduced: value for coefficient A_reduced.
1823 @type A_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1824 @keyword B: value for coefficient B
1825 @type B: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1826 @keyword B_reduced: value for coefficient B_reduced
1827 @type B_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1828 @keyword C: value for coefficient C
1829 @type C: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1830 @keyword C_reduced: value for coefficient C_reduced
1831 @type C_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1832 @keyword D: value for coefficient D
1833 @type D: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1834 @keyword D_reduced: value for coefficient D_reduced
1835 @type D_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1836 @keyword X: value for coefficient X
1837 @type X: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1838 @keyword X_reduced: value for coefficient X_reduced
1839 @type X_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1840 @keyword Y: value for coefficient Y
1841 @type Y: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1842 @keyword Y_reduced: value for coefficient Y_reduced
1843 @type Y_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.Function>}.
1844 @keyword d: value for coefficient d
1845 @type d: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1846 @keyword d_reduced: value for coefficient d_reduced
1847 @type d_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
1848 @keyword y: value for coefficient y
1849 @type y: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1850 @keyword d_contact: value for coefficient d_contact
1851 @type d_contact: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnContactOne<escript.FunctionOnContactOne>} or L{FunctionOnContactZero<escript.FunctionOnContactZero>}.
1852 @keyword d_contact_reduced: value for coefficient d_contact_reduced
1853 @type d_contact_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>} or L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>}.
1854 @keyword y_contact: value for coefficient y_contact
1855 @type y_contact: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnContactOne<escript.FunctionOnContactOne>} or L{FunctionOnContactZero<escript.FunctionOnContactZero>}.
1856 @keyword y_contact_reduced: value for coefficient y_contact_reduced
1857 @type y_contact_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunctionOnContactOne<escript.FunctionOnContactOne>} or L{ReducedFunctionOnContactZero<escript.FunctionOnContactZero>}.
1858 @keyword r: values prescribed to the solution at the locations of constraints
1859 @type r: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1860 depending of reduced order is used for the solution.
1861 @keyword q: mask for location of constraints
1862 @type q: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1863 depending of reduced order is used for the representation of the equation.
1864 @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1865 """
1866 super(LinearPDE,self).setValue(**coefficients)
1867 # check if the systrem is inhomogeneous:
1868 if len(coefficients)>0 and not self.isUsingLumping():
1869 q=self.getCoefficient("q")
1870 r=self.getCoefficient("r")
1871 if not q.isEmpty() and not r.isEmpty():
1872 if util.Lsup(q*r)>0.:
1873 self.trace("Inhomogeneous constraint detected.")
1874 self.invalidateSystem()
1875
1876
1877 def getResidual(self,u=None):
1878 """
1879 return the residual of u or the current solution if u is not present.
1880
1881 @param u: argument in the residual calculation. It must be representable in C{elf.getFunctionSpaceForSolution()}. If u is not present or equals L{None}
1882 the current solution is used.
1883 @type u: L{Data<escript.Data>} or None
1884 @return: residual of u
1885 @rtype: L{Data<escript.Data>}
1886 """
1887 if u==None:
1888 return self.getOperator()*self.getSolution()-self.getRightHandSide()
1889 else:
1890 return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())-self.getRightHandSide()
1891
1892 def getFlux(self,u=None):
1893 """
1894 returns the flux M{J} for a given M{u}
1895
1896 M{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]}
1897
1898 or
1899
1900 M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[l]+(B[j]+B_reduced[j])u-X[j]-X_reduced[j]}
1901
1902 @param u: argument in the flux. If u is not present or equals L{None} the current solution is used.
1903 @type u: L{Data<escript.Data>} or None
1904 @return: flux
1905 @rtype: L{Data<escript.Data>}
1906 """
1907 if u==None: u=self.getSolution()
1908 return util.tensormult(self.getCoefficient("A"),util.grad(u,Funtion(self.getDomain))) \
1909 +util.matrixmult(self.getCoefficient("B"),u) \
1910 -util.self.getCoefficient("X") \
1911 +util.tensormult(self.getCoefficient("A_reduced"),util.grad(u,ReducedFuntion(self.getDomain))) \
1912 +util.matrixmult(self.getCoefficient("B_reduced"),u) \
1913 -util.self.getCoefficient("X_reduced")
1914
1915
1916 class Poisson(LinearPDE):
1917 """
1918 Class to define a Poisson equation problem, which is genear L{LinearPDE} of the form
1919
1920 M{-grad(grad(u)[j])[j] = f}
1921
1922 with natural boundary conditons
1923
1924 M{n[j]*grad(u)[j] = 0 }
1925
1926 and constraints:
1927
1928 M{u=0} where M{q>0}
1929
1930 """
1931
1932 def __init__(self,domain,debug=False):
1933 """
1934 initializes a new Poisson equation
1935
1936 @param domain: domain of the PDE
1937 @type domain: L{Domain<escript.Domain>}
1938 @param debug: if True debug informations are printed.
1939
1940 """
1941 super(Poisson, self).__init__(domain,1,1,debug)
1942 self.introduceCoefficients(
1943 f=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1944 f_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE))
1945 self.setSymmetryOn()
1946
1947 def setValue(self,**coefficients):
1948 """
1949 sets new values to coefficients
1950
1951 @param coefficients: new values assigned to coefficients
1952 @keyword f: value for right hand side M{f}
1953 @type f: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1954 @keyword q: mask for location of constraints
1955 @type q: any type that can be casted to rank zeo L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1956 depending of reduced order is used for the representation of the equation.
1957 @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1958 """
1959 super(Poisson, self).setValue(**coefficients)
1960
1961
1962 def getCoefficient(self,name):
1963 """
1964 return the value of the coefficient name of the general PDE
1965 @param name: name of the coefficient requested.
1966 @type name: C{string}
1967 @return: the value of the coefficient name
1968 @rtype: L{Data<escript.Data>}
1969 @raise IllegalCoefficient: invalid coefficient name
1970 @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
1971 """
1972 if name == "A" :
1973 return escript.Data(util.kronecker(self.getDim()),escript.Function(self.getDomain()))
1974 elif name == "Y" :
1975 return self.getCoefficient("f")
1976 elif name == "Y_reduced" :
1977 return self.getCoefficient("f_reduced")
1978 else:
1979 return super(Poisson, self).getCoefficient(name)
1980
1981 class Helmholtz(LinearPDE):
1982 """
1983 Class to define a Helmhotz equation problem, which is genear L{LinearPDE} of the form
1984
1985 M{S{omega}*u - grad(k*grad(u)[j])[j] = f}
1986
1987 with natural boundary conditons
1988
1989 M{k*n[j]*grad(u)[j] = g- S{alpha}u }
1990
1991 and constraints:
1992
1993 M{u=r} where M{q>0}
1994
1995 """
1996
1997 def __init__(self,domain,debug=False):
1998 """
1999 initializes a new Poisson equation
2000
2001 @param domain: domain of the PDE
2002 @type domain: L{Domain<escript.Domain>}
2003 @param debug: if True debug informations are printed.
2004
2005 """
2006 super(Helmholtz, self).__init__(domain,1,1,debug)
2007 self.introduceCoefficients(
2008 omega=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.OPERATOR),
2009 k=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.OPERATOR),
2010 f=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2011 f_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2012 alpha=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.OPERATOR),
2013 g=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2014 g_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE))
2015 self.setSymmetryOn()
2016
2017 def setValue(self,**coefficients):
2018 """
2019 sets new values to coefficients
2020
2021 @param coefficients: new values assigned to coefficients
2022 @keyword omega: value for coefficient M{S{omega}}
2023 @type omega: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
2024 @keyword k: value for coefficeint M{k}
2025 @type k: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
2026 @keyword f: value for right hand side M{f}
2027 @type f: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
2028 @keyword alpha: value for right hand side M{S{alpha}}
2029 @type alpha: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
2030 @keyword g: value for right hand side M{g}
2031 @type g: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
2032 @keyword r: prescribed values M{r} for the solution in constraints.
2033 @type r: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2034 depending of reduced order is used for the representation of the equation.
2035 @keyword q: mask for location of constraints
2036 @type q: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2037 depending of reduced order is used for the representation of the equation.
2038 @raise IllegalCoefficient: if an unknown coefficient keyword is used.
2039 """
2040 super(Helmholtz, self).setValue(**coefficients)
2041
2042 def getCoefficient(self,name):
2043 """
2044 return the value of the coefficient name of the general PDE
2045
2046 @param name: name of the coefficient requested.
2047 @type name: C{string}
2048 @return: the value of the coefficient name
2049 @rtype: L{Data<escript.Data>}
2050 @raise IllegalCoefficient: invalid name
2051 """
2052 if name == "A" :
2053 return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))*self.getCoefficient("k")
2054 elif name == "D" :
2055 return self.getCoefficient("omega")
2056 elif name == "Y" :
2057 return self.getCoefficient("f")
2058 elif name == "d" :
2059 return self.getCoefficient("alpha")
2060 elif name == "y" :
2061 return self.getCoefficient("g")
2062 elif name == "Y_reduced" :
2063 return self.getCoefficient("f_reduced")
2064 elif name == "y_reduced" :
2065 return self.getCoefficient("g_reduced")
2066 else:
2067 return super(Helmholtz, self).getCoefficient(name)
2068
2069 class LameEquation(LinearPDE):
2070 """
2071 Class to define a Lame equation problem:
2072
2073 M{-grad(S{mu}*(grad(u[i])[j]+grad(u[j])[i]))[j] - grad(S{lambda}*grad(u[k])[k])[j] = F_i -grad(S{sigma}[ij])[j] }
2074
2075 with natural boundary conditons:
2076
2077 M{n[j]*(S{mu}*(grad(u[i])[j]+grad(u[j])[i]) + S{lambda}*grad(u[k])[k]) = f_i +n[j]*S{sigma}[ij] }
2078
2079 and constraints:
2080
2081 M{u[i]=r[i]} where M{q[i]>0}
2082
2083 """
2084
2085 def __init__(self,domain,debug=False):
2086 super(LameEquation, self).__init__(domain,\
2087 domain.getDim(),domain.getDim(),debug)
2088 self.introduceCoefficients(lame_lambda=PDECoef(PDECoef.INTERIOR,(),PDECoef.OPERATOR),
2089 lame_mu=PDECoef(PDECoef.INTERIOR,(),PDECoef.OPERATOR),
2090 F=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2091 sigma=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
2092 f=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE))
2093 self.setSymmetryOn()
2094
2095 def setValues(self,**coefficients):
2096 """
2097 sets new values to coefficients
2098
2099 @param coefficients: new values assigned to coefficients
2100 @keyword lame_mu: value for coefficient M{S{mu}}
2101 @type lame_mu: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
2102 @keyword lame_lambda: value for coefficient M{S{lambda}}
2103 @type lame_lambda: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
2104 @keyword F: value for internal force M{F}
2105 @type F: any type that can be casted to L{Vector<escript.Vector>} object on L{Function<escript.Function>}
2106 @keyword sigma: value for initial stress M{S{sigma}}
2107 @type sigma: any type that can be casted to L{Tensor<escript.Tensor>} object on L{Function<escript.Function>}
2108 @keyword f: value for extrenal force M{f}
2109 @type f: any type that can be casted to L{Vector<escript.Vector>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}
2110 @keyword r: prescribed values M{r} for the solution in constraints.
2111 @type r: any type that can be casted to L{Vector<escript.Vector>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2112 depending of reduced order is used for the representation of the equation.
2113 @keyword q: mask for location of constraints
2114 @type q: any type that can be casted to L{Vector<escript.Vector>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2115 depending of reduced order is used for the representation of the equation.
2116 @raise IllegalCoefficient: if an unknown coefficient keyword is used.
2117 """
2118 super(LameEquation, self).setValues(**coefficients)
2119
2120 def getCoefficient(self,name):
2121 """
2122 return the value of the coefficient name of the general PDE
2123
2124 @param name: name of the coefficient requested.
2125 @type name: C{string}
2126 @return: the value of the coefficient name
2127 @rtype: L{Data<escript.Data>}
2128 @raise IllegalCoefficient: invalid coefficient name
2129 """
2130 if name == "A" :
2131 out =self.createCoefficient("A")
2132 for i in range(self.getDim()):
2133 for j in range(self.getDim()):
2134 out[i,i,j,j] += self.getCoefficient("lame_lambda")
2135 out[i,j,j,i] += self.getCoefficient("lame_mu")
2136 out[i,j,i,j] += self.getCoefficient("lame_mu")
2137 return out
2138 elif name == "X" :
2139 return self.getCoefficient("sigma")
2140 elif name == "Y" :
2141 return self.getCoefficient("F")
2142 elif name == "y" :
2143 return self.getCoefficient("f")
2144 else:
2145 return super(LameEquation, self).getCoefficient(name)
2146
2147 def LinearSinglePDE(domain,debug=False):
2148 """
2149 defines a single linear PDEs
2150
2151 @param domain: domain of the PDE
2152 @type domain: L{Domain<escript.Domain>}
2153 @param debug: if True debug informations are printed.
2154 @rtype: L{LinearPDE}
2155 """
2156 return LinearPDE(domain,numEquations=1,numSolutions=1,debug=debug)
2157
2158 def LinearPDESystem(domain,debug=False):
2159 """
2160 defines a system of linear PDEs
2161
2162 @param domain: domain of the PDE
2163 @type domain: L{Domain<escript.Domain>}
2164 @param debug: if True debug informations are printed.
2165 @rtype: L{LinearPDE}
2166 """
2167 return LinearPDE(domain,numEquations=domain.getDim(),numSolutions=domain.getDim(),debug=debug)
2168
2169 class TransportPDE(LinearProblem):
2170 """
2171 This class is used to define a transport problem given by a general linear, time dependend, second order PDE
2172 for an unknown, non-negative, time-dependend function M{u} on a given domain defined through a L{Domain<escript.Domain>} object.
2173
2174 For a single equation with a solution with a single component the transport problem is defined in the following form:
2175
2176 M{(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)}
2177
2178
2179 where M{u_t} denotes the time derivative of M{u} and M{grad(F)} denotes the spatial derivative of M{F}. Einstein's summation convention,
2180 ie. summation over indexes appearing twice in a term of a sum is performed, is used.
2181 The coefficients M{M}, M{A}, M{B}, M{C}, M{D}, M{X} and M{Y} have to be specified through L{Data<escript.Data>} objects in the
2182 L{Function<escript.Function>} and the coefficients M{M_reduced}, M{A_reduced}, M{B_reduced}, M{C_reduced}, M{D_reduced}, M{X_reduced} and M{Y_reduced} have to be specified through L{Data<escript.Data>} objects in the
2183 L{ReducedFunction<escript.ReducedFunction>}. It is also allowed to use objects that can be converted into
2184 such L{Data<escript.Data>} objects. M{M} and M{M_reduced} are scalar, M{A} and M{A_reduced} are rank two,
2185 M{B}, M{C}, M{X}, M{B_reduced}, M{C_reduced} and M{X_reduced} are rank one and M{D}, M{D_reduced}, M{Y} and M{Y_reduced} are scalar.
2186
2187 The following natural boundary conditions are considered:
2188
2189 M{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}
2190
2191 where M{n} is the outer normal field. Notice that the coefficients M{A}, M{A_reduced}, M{B}, M{B_reduced}, M{X} and M{X_reduced} are defined in the transport problem. The coefficients M{m}, M{d} and M{y} and are each a scalar in the L{FunctionOnBoundary<escript.FunctionOnBoundary>} and the coefficients M{m_reduced}, M{d_reduced} and M{y_reduced} and are each a scalar in the L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
2192
2193 Constraints for the solution prescribing the value of the solution at certain locations in the domain. They have the form
2194
2195 M{u_t=r} where M{q>0}
2196
2197 M{r} and M{q} are each scalar where M{q} is the characteristic function defining where the constraint is applied.
2198 The constraints override any other condition set by the transport problem or the boundary condition.
2199
2200 The transport problem is symmetrical if
2201
2202 M{A[i,j]=A[j,i]} and M{B[j]=C[j]} and M{A_reduced[i,j]=A_reduced[j,i]} and M{B_reduced[j]=C_reduced[j]}
2203
2204 For a system and a solution with several components the transport problem has the form
2205
2206 M{(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] }
2207
2208 M{A} and M{A_reduced} are of rank four, M{B}, M{B_reduced}, M{C} and M{C_reduced} are each of rank three, M{M}, M{M_reduced}, M{D}, M{D_reduced}, M{X_reduced} and M{X} are each a rank two and M{Y} and M{Y_reduced} are of rank one. The natural boundary conditions take the form:
2209
2210 M{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}
2211
2212 The coefficient M{d} and M{m} are of rank two and M{y} are of rank one with all in the L{FunctionOnBoundary<escript.FunctionOnBoundary>}. Constraints take the form and the coefficients M{d_reduced} and M{m_reduced} are of rank two and M{y_reduced} is of rank one with all in the L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
2213
2214 Constraints take the form
2215
2216 M{u[i]_t=r[i]} where M{q[i]>0}
2217
2218 M{r} and M{q} are each rank one. Notice that at some locations not necessarily all components must have a constraint.
2219
2220 The transport problem is symmetrical if
2221
2222 - M{M[i,k]=M[i,k]}
2223 - M{M_reduced[i,k]=M_reduced[i,k]}
2224 - M{A[i,j,k,l]=A[k,l,i,j]}
2225 - M{A_reduced[i,j,k,l]=A_reduced[k,l,i,j]}
2226 - M{B[i,j,k]=C[k,i,j]}
2227 - M{B_reduced[i,j,k]=C_reduced[k,i,j]}
2228 - M{D[i,k]=D[i,k]}
2229 - M{D_reduced[i,k]=D_reduced[i,k]}
2230 - M{m[i,k]=m[k,i]}
2231 - M{m_reduced[i,k]=m_reduced[k,i]}
2232 - M{d[i,k]=d[k,i]}
2233 - M{d_reduced[i,k]=d_reduced[k,i]}
2234
2235 L{TransportPDE} also supports solution discontinuities over a contact region in the domain. To specify the conditions across the
2236 discontinuity we are using the generalised flux M{J} which is in the case of a systems of PDEs and several components of the solution
2237 defined as
2238
2239 M{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]}
2240
2241 For the case of single solution component and single PDE M{J} is defined
2242
2243 M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[j]+(B[i]+B_reduced[i])*u+X[i]+X_reduced[i]}
2244
2245 In the context of discontinuities M{n} denotes the normal on the discontinuity pointing from side 0 towards side 1
2246 calculated from L{getNormal<escript.FunctionSpace.getNormal>} of L{FunctionOnContactZero<escript.FunctionOnContactZero>}. For a system of transport problems the contact condition takes the form
2247
2248 M{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]}
2249
2250 where M{J0} and M{J1} are the fluxes on side 0 and side 1 of the discontinuity, respectively. M{jump(u)}, which is the difference
2251 of the solution at side 1 and at side 0, denotes the jump of M{u} across discontinuity along the normal calculated by
2252 L{jump<util.jump>}.
2253 The coefficient M{d_contact} is a rank two and M{y_contact} is a rank one both in the L{FunctionOnContactZero<escript.FunctionOnContactZero>} or L{FunctionOnContactOne<escript.FunctionOnContactOne>}.
2254 The coefficient M{d_contact_reduced} is a rank two and M{y_contact_reduced} is a rank one both in the L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}.
2255 In case of a single PDE and a single component solution the contact condition takes the form
2256
2257 M{n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)}
2258
2259 In this case the coefficient M{d_contact} and M{y_contact} are each scalar both in the L{FunctionOnContactZero<escript.FunctionOnContactZero>} or L{FunctionOnContactOne<escript.FunctionOnContactOne>} and the coefficient M{d_contact_reduced} and M{y_contact_reduced} are each scalar both in the L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}
2260
2261 typical usage:
2262
2263 p=TransportPDE(dom)
2264 p.setValue(M=1.,C=[-1.,0.])
2265 p.setInitialSolution(u=exp(-length(dom.getX()-[0.1,0.1])**2)
2266 t=0
2267 dt=0.1
2268 while (t<1.):
2269 u=p.solve(dt)
2270
2271 """
2272 def __init__(self,domain,num_equations=None,numSolutions=None, theta=0.5,debug=False):
2273 """
2274 initializes a linear problem
2275
2276 @param domain: domain of the PDE
2277 @type domain: L{Domain<escript.Domain>}
2278 @param numEquations: number of equations. If numEquations==None the number of equations
2279 is extracted from the coefficients.
2280 @param numSolutions: number of solution components. If numSolutions==None the number of solution components
2281 is extracted from the coefficients.
2282 @param debug: if True debug informations are printed.
2283 @param theta: time stepping control: =1 backward Euler, =0 forward Euler, =0.5 Crank-Nicholson scheme, U{see here<http://en.wikipedia.org/wiki/Crank-Nicolson_method>}.
2284
2285 """
2286 if theta<0 or theta>1:
2287 raise ValueError,"theta (=%s) needs to be between 0 and 1 (inclusive)."
2288 super(TransportPDE, self).__init__(domain,num_equations,numSolutions,debug)
2289
2290 self.__theta=theta
2291 #
2292 # the coefficients of the transport problem
2293 #
2294 self.introduceCoefficients(
2295 M=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2296 A=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2297 B=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2298 C=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2299 D=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2300 X=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
2301 Y=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2302 m=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2303 d=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2304 y=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2305 d_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2306 y_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2307 M_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2308 A_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2309 B_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2310 C_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2311 D_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2312 X_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
2313 Y_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2314 m_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2315 d_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2316 y_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2317 d_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2318 y_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2319 r=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.RIGHTHANDSIDE),
2320 q=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.BOTH) )
2321
2322 def __str__(self):
2323 """
2324 returns string representation of the transport problem
2325
2326 @return: a simple representation of the transport problem
2327 @rtype: C{str}
2328 """
2329 return "<TransportPDE %d>"%id(self)
2330
2331 def getTheta(self):
2332 """
2333 returns the time stepping control parameter
2334 @rtype: float
2335 """
2336 return self.__theta
2337 def getRequiredSystemType(self):
2338 """
2339 returns the system type which needs to be used by the current set up.
2340 """
2341 return 0
2342 return self.getDomain().getTransportMatrixTypeId(self.getSolverMethod()[0],self.getSolverMethod()[1],self.getSolverPackage(),self.isSymmetric())
2343
2344 def checkSymmetry(self,verbose=True):
2345 """
2346 test the transport problem for symmetry.
2347
2348 @param verbose: if equal to True or not present a report on coefficients which are breaking the symmetry is printed.
2349 @type verbose: C{bool}
2350 @return: True if the PDE is symmetric.
2351 @rtype: C{bool}
2352 @note: This is a very expensive operation. It should be used for degugging only! The symmetry flag is not altered.
2353 """
2354 out=True
2355 out=out and self.checkSymmetricTensor("M", verbose)
2356 out=out and self.checkSymmetricTensor("M_reduced", verbose)
2357 out=out and self.checkSymmetricTensor("A", verbose)
2358 out=out and self.checkSymmetricTensor("A_reduced", verbose)
2359 out=out and self.checkReciprocalSymmetry("B","C", verbose)
2360 out=out and self.checkReciprocalSymmetry("B_reduced","C_reduced", verbose)
2361 out=out and self.checkSymmetricTensor("D", verbose)
2362 out=out and self.checkSymmetricTensor("D_reduced", verbose)
2363 out=out and self.checkSymmetricTensor("m", verbose)
2364 out=out and self.checkSymmetricTensor("m_reduced", verbose)
2365 out=out and self.checkSymmetricTensor("d", verbose)
2366 out=out and self.checkSymmetricTensor("d_reduced", verbose)
2367 out=out and self.checkSymmetricTensor("d_contact", verbose)
2368 out=out and self.checkSymmetricTensor("d_contact_reduced", verbose)
2369 return out
2370
2371 def setValue(self,**coefficients):
2372 """
2373 sets new values to coefficients
2374
2375 @param coefficients: new values assigned to coefficients
2376 @keyword M: value for coefficient M.
2377 @type M: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
2378 @keyword M_reduced: value for coefficient M_reduced.
2379 @type M_reduced: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.ReducedFunction>}.
2380 @keyword A: value for coefficient A.
2381 @type A: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
2382 @keyword A_reduced: value for coefficient A_reduced.
2383 @type A_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
2384 @keyword B: value for coefficient B
2385 @type B: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
2386 @keyword B_reduced: value for coefficient B_reduced
2387 @type B_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
2388 @keyword C: value for coefficient C
2389 @type C: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
2390 @keyword C_reduced: value for coefficient C_reduced
2391 @type C_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
2392 @keyword D: value for coefficient D
2393 @type D: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
2394 @keyword D_reduced: value for coefficient D_reduced
2395 @type D_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
2396 @keyword X: value for coefficient X
2397 @type X: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
2398 @keyword X_reduced: value for coefficient X_reduced
2399 @type X_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
2400 @keyword Y: value for coefficient Y
2401 @type Y: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
2402 @keyword Y_reduced: value for coefficient Y_reduced
2403 @type Y_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.Function>}.
2404 @keyword m: value for coefficient m
2405 @type m: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
2406 @keyword m_reduced: value for coefficient m_reduced
2407 @type m_reduced: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
2408 @keyword d: value for coefficient d
2409 @type d: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
2410 @keyword d_reduced: value for coefficient d_reduced
2411 @type d_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
2412 @keyword y: value for coefficient y
2413 @type y: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
2414 @keyword d_contact: value for coefficient d_contact
2415 @type d_contact: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnContactOne<escript.FunctionOnContactOne>} or L{FunctionOnContactZero<escript.FunctionOnContactZero>}.
2416 @keyword d_contact_reduced: value for coefficient d_contact_reduced
2417 @type d_contact_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>} or L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>}.
2418 @keyword y_contact: value for coefficient y_contact
2419 @type y_contact: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnContactOne<escript.FunctionOnContactOne>} or L{FunctionOnContactZero<escript.FunctionOnContactZero>}.
2420 @keyword y_contact_reduced: value for coefficient y_contact_reduced
2421 @type y_contact_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunctionOnContactOne<escript.FunctionOnContactOne>} or L{ReducedFunctionOnContactZero<escript.FunctionOnContactZero>}.
2422 @keyword r: values prescribed to the solution at the locations of constraints
2423 @type r: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2424 depending of reduced order is used for the solution.
2425 @keyword q: mask for location of constraints
2426 @type q: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2427 depending of reduced order is used for the representation of the equation.
2428 @raise IllegalCoefficient: if an unknown coefficient keyword is used.
2429 """
2430 super(TransportPDE,self).setValue(**coefficients)
2431
2432 def createOperator(self):
2433 """
2434 returns an instance of a new transport operator
2435 """
2436 self.trace("New Transport problem is allocated.")
2437 return self.getDomain().newTransportProblem( \
2438 self.getTheta(),
2439 self.getNumEquations(), \
2440 self.getFunctionSpaceForSolution(), \
2441 self.getSystemType())
2442
2443 def getSolution(self,dt,**options):
2444 """
2445 returns the solution of the problem
2446 @return: the solution
2447 @rtype: L{Data<escript.Data>}
2448
2449 @note: This method is overwritten when implementing a particular linear problem.
2450 """
2451 if not self.isSolutionValid():
2452 options[self.TOLERANCE_KEY]=self.getTolerance()
2453 self.setSolution(self.getOperator().solve(self.getRightHandSide(),dt,options))
2454 self.validSolution()
2455 return self.getCurrentSolution()
2456 def setInitialSolution(self,u):
2457 self.getOperator().setInitialValue(util.interpolate(u,self.getFunctionSpaceForSolution()))
2458 def getSafeTimeStepSize(self):
2459 return self.getOperator().getSafeTimeStepSize()
2460
2461 def getSystem(self):
2462 """
2463 return the operator and right hand side of the PDE
2464
2465 @return: the discrete version of the PDE
2466 @rtype: C{tuple} of L{Operator,<escript.Operator>} and L{Data<escript.Data>}.
2467
2468 @note: This method is overwritten when implementing a particular linear problem.
2469 """
2470 if not self.isOperatorValid() or not self.isRightHandSideValid():
2471 self.resetRightHandSide()
2472 righthandside=self.getCurrentRightHandSide()
2473 self.resetOperator()
2474 operator=self.getCurrentOperator()
2475 self.getDomain().addPDEToTransportProblem(
2476 operator,
2477 righthandside,
2478 self.getCoefficient("M"),
2479 self.getCoefficient("A"),
2480 self.getCoefficient("B"),
2481 self.getCoefficient("C"),
2482 self.getCoefficient("D"),
2483 self.getCoefficient("X"),
2484 self.getCoefficient("Y"),
2485 self.getCoefficient("d"),
2486 self.getCoefficient("y"),
2487 self.getCoefficient("d_contact"),
2488 self.getCoefficient("y_contact"))
2489 operator.insertConstraint(righthandside,self.getCoefficient("q"),self.getCoefficient("r"))
2490 self.trace("New system has been built.")
2491 self.validOperator()
2492 self.validRightHandSide()
2493 return (self.getCurrentOperator(), self.getCurrentRightHandSide())

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26