/[escript]/trunk-mpi-branch/escript/py_src/linearPDEs.py
ViewVC logotype

Contents of /trunk-mpi-branch/escript/py_src/linearPDEs.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 148 - (show annotations)
Tue Aug 23 01:24:31 2005 UTC (14 years, 1 month ago) by jgs
Original Path: trunk/esys2/escript/py_src/linearPDEs.py
File MIME type: text/x-python
File size: 68376 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-08-23

1 # $Id$
2
3 ## @file linearPDEs.py
4
5 """
6 Functions and classes for linear PDEs
7 """
8
9 import escript
10 import util
11 import numarray
12
13
14 class IllegalCoefficient(ValueError):
15 """
16 raised if an illegal coefficient of the general ar particular PDE is requested.
17 """
18
19 class IllegalCoefficientValue(ValueError):
20 """
21 raised if an incorrect value for a coefficient is used.
22 """
23
24 class UndefinedPDEError(ValueError):
25 """
26 raised if a PDE is not fully defined yet.
27 """
28
29 def _CompTuple2(t1,t2):
30 """
31 Compare two tuples
32
33 @param t1 The first tuple
34 @param t2 The second tuple
35
36 """
37
38 dif=t1[0]+t1[1]-(t2[0]+t2[1])
39 if dif<0: return 1
40 elif dif>0: return -1
41 else: return 0
42
43 class PDECoefficient:
44 """
45 A class for PDE coefficients
46 """
47 # identifier for location of Data objects defining COEFFICIENTS
48 INTERIOR=0
49 BOUNDARY=1
50 CONTACT=2
51 CONTINUOUS=3
52 # identifier in the pattern of COEFFICIENTS:
53 # the pattern is a tuple of EQUATION,SOLUTION,DIM where DIM represents the spatial dimension, EQUATION the number of equations and SOLUTION the
54 # number of unknowns.
55 EQUATION=3
56 SOLUTION=4
57 DIM=5
58 # indicator for what is altered if the coefficient is altered:
59 OPERATOR=5
60 RIGHTHANDSIDE=6
61 BOTH=7
62 def __init__(self,where,pattern,altering):
63 """
64 Initialise a PDE Coefficient type
65 """
66 self.what=where
67 self.pattern=pattern
68 self.altering=altering
69 self.resetValue()
70
71 def resetValue(self):
72 """
73 resets coefficient value to default
74 """
75 self.value=escript.Data()
76
77 def getFunctionSpace(self,domain):
78 """
79 defines the FunctionSpace of the coefficient on the domain
80
81 @param domain:
82 """
83 if self.what==self.INTERIOR: return escript.Function(domain)
84 elif self.what==self.BOUNDARY: return escript.FunctionOnBoundary(domain)
85 elif self.what==self.CONTACT: return escript.FunctionOnContactZero(domain)
86 elif self.what==self.CONTINUOUS: return escript.ContinuousFunction(domain)
87
88 def getValue(self):
89 """
90 returns the value of the coefficient:
91 """
92 return self.value
93
94 def setValue(self,domain,numEquations=1,numSolutions=1,newValue=None):
95 """
96 set the value of the coefficient to new value
97 """
98 if newValue==None:
99 newValue=escript.Data()
100 elif isinstance(newValue,escript.Data):
101 if not newValue.isEmpty():
102 newValue=escript.Data(newValue,self.getFunctionSpace(domain))
103 else:
104 newValue=escript.Data(newValue,self.getFunctionSpace(domain))
105 if not newValue.isEmpty():
106 if not self.getShape(domain,numEquations,numSolutions)==newValue.getShape():
107 raise IllegalCoefficientValue,"Expected shape for coefficient %s is %s but actual shape is %s."%(self.getShape(domain,numEquations,numSolutions),newValue.getShape())
108 self.value=newValue
109
110 def isAlteringOperator(self):
111 """
112 return true if the operator of the PDE is changed when the coefficient is changed
113 """
114 if self.altering==self.OPERATOR or self.altering==self.BOTH:
115 return not None
116 else:
117 return None
118
119 def isAlteringRightHandSide(self):
120 """
121 return true if the right hand side of the PDE is changed when the coefficient is changed
122 """
123 if self.altering==self.RIGHTHANDSIDE or self.altering==self.BOTH:
124 return not None
125 else:
126 return None
127
128 def estimateNumEquationsAndNumSolutions(self,domain,shape=()):
129 """
130 tries to estimate the number of equations in a given tensor shape for a given spatial dimension dim
131
132 @param shape:
133 @param dim:
134 """
135 dim=domain.getDim()
136 if len(shape)>0:
137 num=max(shape)+1
138 else:
139 num=1
140 search=[]
141 for u in range(num):
142 for e in range(num):
143 search.append((e,u))
144 search.sort(_CompTuple2)
145 for item in search:
146 s=self.getShape(domain,item[0],item[1])
147 if len(s)==0 and len(shape)==0:
148 return (1,1)
149 else:
150 if s==shape: return item
151 return None
152
153 def getShape(self,domain,numEquations=1,numSolutions=1):
154 """
155 builds the required shape for a given number of equations e, number of unknowns u and spatial dimension dim
156
157 @param e:
158 @param u:
159 @param dim:
160 """
161 dim=domain.getDim()
162 s=()
163 for i in self.pattern:
164 if i==self.EQUATION:
165 if numEquations>1: s=s+(numEquations,)
166 elif i==self.SOLUTION:
167 if numSolutions>1: s=s+(numSolutions,)
168 else:
169 s=s+(dim,)
170 return s
171
172 class LinearPDE:
173 """
174 Class to define a linear PDE of the form
175
176 \f[
177 -(A_{ijkl}u_{k,l})_{,j} -(B_{ijk}u_k)_{,j} + C_{ikl}u_{k,l} +D_{ik}u_k = - (X_{ij})_{,j} + Y_i
178 \f]
179
180 with boundary conditons:
181
182 \f[
183 n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_{ik}u_k = - n_j*X_{ij} + y_i
184 \f]
185
186 and contact conditions
187
188 \f[
189 n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_contact_{ik}[u_k] = - n_j*X_{ij} + y_contact_i
190 \f]
191
192 and constraints:
193
194 \f[
195 u_i=r_i \quad \mathrm{where} \quad q_i>0
196 \f]
197
198 """
199 TOL=1.e-13
200 # solver methods
201 UNKNOWN=-1
202 DEFAULT_METHOD=0
203 DIRECT=1
204 CHOLEVSKY=2
205 PCG=3
206 CR=4
207 CGS=5
208 BICGSTAB=6
209 SSOR=7
210 ILU0=8
211 ILUT=9
212 JACOBI=10
213 GMRES=11
214 PRES20=12
215 LUMPING=13
216 # matrix reordering:
217 NO_REORDERING=0
218 MINIMUM_FILL_IN=1
219 NESTED_DISSECTION=2
220 # important keys in the dictonary used to hand over options to the solver library:
221 METHOD_KEY="method"
222 SYMMETRY_KEY="symmetric"
223 TOLERANCE_KEY="tolerance"
224
225
226 def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):
227 """
228 initializes a new linear PDE
229
230 @param domain: domain of the PDE
231 @type domain: L{Domain}
232 @param numEquations: number of equations. If numEquations==None the number of equations
233 is exracted from the PDE coefficients.
234 @param numSolutions: number of solution components. If numSolutions==None the number of solution components
235 is exracted from the PDE coefficients.
236 @param debug: if True debug informations are printed.
237
238
239 """
240 #
241 # the coefficients of the general PDE:
242 #
243 self.__COEFFICIENTS_OF_GENEARL_PDE={
244 "A" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM,PDECoefficient.SOLUTION,PDECoefficient.DIM),PDECoefficient.OPERATOR),
245 "B" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),
246 "C" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION,PDECoefficient.DIM),PDECoefficient.OPERATOR),
247 "D" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),
248 "X" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM),PDECoefficient.RIGHTHANDSIDE),
249 "Y" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
250 "d" : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),
251 "y" : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
252 "d_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),
253 "y_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
254 "r" : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
255 "q" : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.SOLUTION,),PDECoefficient.BOTH)}
256
257 # COEFFICIENTS can be overwritten by subclasses:
258 self.COEFFICIENTS=self.__COEFFICIENTS_OF_GENEARL_PDE
259 # initialize attributes
260 self.__debug=debug
261 self.__domain=domain
262 self.__numEquations=numEquations
263 self.__numSolutions=numSolutions
264 self.__resetSystem()
265
266 # set some default values:
267 self.__homogeneous_constraint=True
268 self.__row_function_space=escript.Solution(self.__domain)
269 self.__column_function_space=escript.Solution(self.__domain)
270 self.__tolerance=1.e-8
271 self.__solver_method=self.DEFAULT_METHOD
272 self.__matrix_type=self.__domain.getSystemMatrixTypeId(self.DEFAULT_METHOD,False)
273 self.__sym=False
274
275 self.resetCoefficients()
276 self.trace("PDE Coeffients are %s"%str(self.COEFFICIENTS.keys()))
277 # =============================================================================
278 # general stuff:
279 # =============================================================================
280 def __str__(self):
281 return "<LinearPDE %d>"%id(self)
282 # =============================================================================
283 # debug :
284 # =============================================================================
285 def setDebugOn(self):
286 """
287 switches on debugging
288 """
289 self.__debug=not None
290
291 def setDebugOff(self):
292 """
293 switches off debugging
294 """
295 self.__debug=None
296
297 def trace(self,text):
298 """
299 print the text message if debugging is swiched on.
300
301 @param name: name of the coefficient enquired.
302 @type name: C{string}
303 """
304 if self.__debug: print "%s: %s"%(str(self),text)
305
306 # =============================================================================
307 # some service functions:
308 # =============================================================================
309 def getDomain(self):
310 """
311 returns the domain of the PDE
312
313 @return : the domain of the PDE
314 @rtype : C{Domain}
315
316 """
317 return self.__domain
318
319 def getDim(self):
320 """
321 returns the spatial dimension of the PDE
322
323 @return : the spatial dimension of the PDE domain
324 @rtype : C{int}
325 """
326 return self.getDomain().getDim()
327
328 def getNumEquations(self):
329 """
330 returns the number of equations
331
332 @return : the number of equations
333 @rtype : C{int}
334 @raise UndefinedPDEError: if the number of equations is not be specified yet.
335 """
336 if self.__numEquations==None:
337 raise UndefinedPDEError,"Number of equations is undefined. Please specify argument numEquations."
338 else:
339 return self.__numEquations
340
341 def getNumSolutions(self):
342 """
343 returns the number of unknowns
344
345 @return : the number of unknowns
346 @rtype : C{int}
347 @raise UndefinedPDEError: if the number of unknowns is not be specified yet.
348 """
349 if self.__numSolutions==None:
350 raise UndefinedPDEError,"Number of solution is undefined. Please specify argument numSolutions."
351 else:
352 return self.__numSolutions
353
354 def getFunctionSpaceForEquation(self):
355 """
356 returns the L{escript.FunctionSpace} used to discretize the equation
357
358 @return : representation space of equation
359 @rtype : L{escript.FunctionSpace}
360
361 """
362 return self.__row_function_space
363
364 def getFunctionSpaceForSolution(self):
365 """
366 returns the L{escript.FunctionSpace} used to represent the solution
367
368 @return : representation space of solution
369 @rtype : L{escript.FunctionSpace}
370
371 """
372 return self.__column_function_space
373
374
375 def getOperator(self):
376 """
377 provides access to the operator of the PDE
378
379 @return : the operator of the PDE
380 @rtype : L{Operator}
381 """
382 m=self.getSystem()[0]
383 if self.isUsingLumping():
384 return self.copyConstraint(1./m)
385 else:
386 return m
387
388 def getRightHandSide(self):
389 """
390 provides access to the right hand side of the PDE
391
392 @return : the right hand side of the PDE
393 @rtype : L{escript.Data}
394 """
395 r=self.getSystem()[1]
396 if self.isUsingLumping():
397 return self.copyConstraint(r)
398 else:
399 return r
400
401 def applyOperator(self,u=None):
402 """
403 applies the operator of the PDE to a given u or the solution of PDE if u is not present.
404
405 @param u: argument of the operator. It must be representable in C{elf.getFunctionSpaceForSolution()}. If u is not present or equals L{None}
406 the current solution is used.
407 @type u: L{escript.Data} or None
408 @return : image of u
409 @rtype : L{escript.Data}
410 """
411 if u==None:
412 return self.getOperator()*self.getSolution()
413 else:
414 self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())
415
416 def getResidual(self,u=None):
417 """
418 return the residual of u or the current solution if u is not present.
419
420 @param u: argument in the residual calculation. It must be representable in C{elf.getFunctionSpaceForSolution()}. If u is not present or equals L{None}
421 the current solution is used.
422 @type u: L{escript.Data} or None
423 @return : residual of u
424 @rtype : L{escript.Data}
425 """
426 return self.applyOperator(u)-self.getRightHandSide()
427
428 def checkSymmetry(self,verbose=True):
429 """
430 test the PDE for symmetry.
431
432
433 @param verbose: if equal to True or not present a report on coefficients which are breaking the symmetry is printed.
434 @type verbose: C{bool}
435 @return: True if the PDE is symmetric.
436 @rtype : C{escript.Data}
437
438 @note: This is a very expensive operation. It should be used for degugging only! The symmetry flag is not altered.
439 """
440 verbose=verbose or self.debug()
441 out=True
442 if self.getNumSolutions()!=self.getNumEquations():
443 if verbose: print "non-symmetric PDE because of different number of equations and solutions"
444 out=False
445 else:
446 A=self.getCoefficientOfGeneralPDE("A")
447 if not A.isEmpty():
448 tol=util.Lsup(A)*self.TOL
449 if self.getNumSolutions()>1:
450 for i in range(self.getNumEquations()):
451 for j in range(self.getDim()):
452 for k in range(self.getNumSolutions()):
453 for l in range(self.getDim()):
454 if util.Lsup(A[i,j,k,l]-A[k,l,i,j])>tol:
455 if verbose: print "non-symmetric PDE because A[%d,%d,%d,%d]!=A[%d,%d,%d,%d]"%(i,j,k,l,k,l,i,j)
456 out=False
457 else:
458 for j in range(self.getDim()):
459 for l in range(self.getDim()):
460 if util.Lsup(A[j,l]-A[l,j])>tol:
461 if verbose: print "non-symmetric PDE because A[%d,%d]!=A[%d,%d]"%(j,l,l,j)
462 out=False
463 B=self.getCoefficientOfGeneralPDE("B")
464 C=self.getCoefficientOfGeneralPDE("C")
465 if B.isEmpty() and not C.isEmpty():
466 if verbose: print "non-symmetric PDE because B is not present but C is"
467 out=False
468 elif not B.isEmpty() and C.isEmpty():
469 if verbose: print "non-symmetric PDE because C is not present but B is"
470 out=False
471 elif not B.isEmpty() and not C.isEmpty():
472 tol=(util.Lsup(B)+util.Lsup(C))*self.TOL/2.
473 if self.getNumSolutions()>1:
474 for i in range(self.getNumEquations()):
475 for j in range(self.getDim()):
476 for k in range(self.getNumSolutions()):
477 if util.Lsup(B[i,j,k]-C[k,i,j])>tol:
478 if verbose: print "non-symmetric PDE because B[%d,%d,%d]!=C[%d,%d,%d]"%(i,j,k,k,i,j)
479 out=False
480 else:
481 for j in range(self.getDim()):
482 if util.Lsup(B[j]-C[j])>tol:
483 if verbose: print "non-symmetric PDE because B[%d]!=C[%d]"%(j,j)
484 out=False
485 if self.getNumSolutions()>1:
486 D=self.getCoefficientOfGeneralPDE("D")
487 if not D.isEmpty():
488 tol=util.Lsup(D)*self.TOL
489 for i in range(self.getNumEquations()):
490 for k in range(self.getNumSolutions()):
491 if util.Lsup(D[i,k]-D[k,i])>tol:
492 if verbose: print "non-symmetric PDE because D[%d,%d]!=D[%d,%d]"%(i,k,k,i)
493 out=False
494
495 return out
496
497 def getSolution(self,**options):
498 """
499 returns the solution of the PDE. If the solution is not valid the PDE is solved.
500
501 @return: the solution
502 @rtype: L{escript.Data}
503 @param options: solver options
504 @keyword verbose:
505 @keyword reordering: reordering scheme to be used during elimination
506 @keyword preconditioner: preconditioner method to be used
507 @keyword iter_max: maximum number of iteration steps allowed.
508 @keyword drop_tolerance:
509 @keyword drop_storage:
510 @keyword truncation:
511 @keyword restart:
512 """
513 if not self.__solution_isValid:
514 mat,f=self.getSystem()
515 if self.isUsingLumping():
516 self.__solution=self.copyConstraint(f*mat)
517 else:
518 options[self.TOLERANCE_KEY]=self.getTolerance()
519 options[self.METHOD_KEY]=self.getSolverMethod()
520 options[self.SYMMETRY_KEY]=self.isSymmetric()
521 self.trace("PDE is resolved.")
522 self.trace("solver options: %s"%str(options))
523 self.__solution=mat.solve(f,options)
524 self.__solution_isValid=True
525 return self.__solution
526
527 def getFlux(self,u=None):
528 """
529 returns the flux J_ij for a given u
530
531 \f[
532 J_ij=A_{ijkl}u_{k,l}+B_{ijk}u_k-X_{ij}
533 \f]
534
535 @param u: argument in the flux. If u is not present or equals L{None} the current solution is used.
536 @type u: L{escript.Data} or None
537 @return : flux
538 @rtype : L{escript.Data}
539
540 """
541 if u==None: u=self.getSolution()
542 return util.tensormult(self.getCoefficientOfGeneralPDE("A"),util.grad(u))+util.matrixmult(self.getCoefficientOfGeneralPDE("B"),u)-util.self.getCoefficientOfGeneralPDE("X")
543
544 # =============================================================================
545 # solver settings:
546 # =============================================================================
547 def setSolverMethod(self,solver=None):
548 """
549 sets a new solver
550
551 @param solver: sets a new solver method.
552 @type solver: C{int}
553
554 """
555 if solver==None: solve=self.DEFAULT_METHOD
556 if not solver==self.getSolverMethod():
557 self.__solver_method=solver
558 self.__checkMatrixType()
559 self.trace("New solver is %s"%self.getSolverMethodName())
560
561 def getSolverMethodName(self):
562 """
563 returns the name of the solver currently used
564
565 @return : the name of the solver currently used.
566 @rtype: C{string}
567 """
568
569 m=self.getSolverMethod()
570 if m==self.DEFAULT_METHOD: return "DEFAULT_METHOD"
571 elif m==self.DIRECT: return "DIRECT"
572 elif m==self.CHOLEVSKY: return "CHOLEVSKY"
573 elif m==self.PCG: return "PCG"
574 elif m==self.CR: return "CR"
575 elif m==self.CGS: return "CGS"
576 elif m==self.BICGSTAB: return "BICGSTAB"
577 elif m==self.SSOR: return "SSOR"
578 elif m==self.GMRES: return "GMRES"
579 elif m==self.PRES20: return "PRES20"
580 elif m==self.LUMPING: return "LUMPING"
581 return None
582
583
584 def getSolverMethod(self):
585 """
586 returns the solver method
587
588 @return : the solver method currently be used.
589 @rtype : C{int}
590 """
591 return self.__solver_method
592
593 def isUsingLumping(self):
594 """
595 checks if matrix lumping is used a solver method
596
597 @return : True is lumping is currently used a solver method.
598 @rtype: C{bool}
599 """
600 return self.getSolverMethod()==self.LUMPING
601
602 def setTolerance(self,tol=1.e-8):
603 """
604 resets the tolerance for the solver method to tol where for an appropriate norm |.|
605
606 |self.getResidual()|<tol*|self.getRightHandSide()|
607
608 defines the stopping criterion.
609
610 @param tol: new tolerance for the solver. If the tol is lower then the current tolerence
611 the system will be resolved.
612 @type solver: C{float}
613 @raise ValueException: if tolerance is not positive.
614 """
615 if not tol>0:
616 raise ValueException,"Tolerance as to be positive"
617 if tol<self.getTolerance(): self.__invalidateSolution()
618 self.trace("New tolerance %e"%tol)
619 self.__tolerance=tol
620 return
621
622 def getTolerance(self):
623 """
624 returns the tolerance set for the solution
625
626 @return: tolerance currently used.
627 @rtype: C{float}
628 """
629 return self.__tolerance
630
631 # =============================================================================
632 # symmetry flag:
633 # =============================================================================
634 def isSymmetric(self):
635 """
636 checks if symmetry is indicated.
637
638 @return : True is a symmetric PDE is indicated, otherwise False is returned
639 @rtype : C{bool}
640 """
641 return self.__sym
642
643 def setSymmetryOn(self):
644 """
645 sets the symmetry flag.
646 """
647 if not self.isSymmetric():
648 self.trace("PDE is set to be symmetric")
649 self.__sym=True
650 self.__checkMatrixType()
651
652 def setSymmetryOff(self):
653 """
654 removes the symmetry flag.
655 """
656 if self.isSymmetric():
657 self.trace("PDE is set to be unsymmetric")
658 self.__sym=False
659 self.__checkMatrixType()
660
661 def setSymmetryTo(self,flag=False):
662 """
663 sets the symmetry flag to flag
664
665 @param flag: If flag, the symmetry flag is set otherwise the symmetry flag is released.
666 @type flag: C{bool}
667 """
668 if flag:
669 self.setSymmetryOn()
670 else:
671 self.setSymmetryOff()
672
673
674 # =============================================================================
675 # function space handling for the equation as well as the solution
676 # =============================================================================
677 def setReducedOrderOn(self):
678 """
679 switches on reduced order for solution and equation representation
680 """
681 self.setReducedOrderForSolutionOn()
682 self.setReducedOrderForEquationOn()
683
684 def setReducedOrderOff(self):
685 """
686 switches off reduced order for solution and equation representation
687 """
688 self.setReducedOrderForSolutionOff()
689 self.setReducedOrderForEquationOff()
690
691 def setReducedOrderTo(self,flag=False):
692 """
693 sets order reduction for both solution and equation representation according to flag.
694
695 @param flag: if flag is True, the order reduction is switched on for both solution and equation representation, otherwise or
696 if flag is not present order reduction is switched off
697 @type flag: C{bool}
698 """
699 self.setReducedOrderForSolutionTo(flag)
700 self.setReducedOrderForEquationTo(flag)
701
702
703 def setReducedOrderForSolutionOn(self):
704 """
705 switches on reduced order for solution representation
706 """
707 new_fs=escript.ReducedSolution(self.getDomain())
708 if self.getFunctionSpaceForSolution()!=new_fs:
709 self.trace("Reduced order is used to solution representation.")
710 self.__column_function_space=new_fs
711 self.__resetSystem()
712
713 def setReducedOrderForSolutionOff(self):
714 """
715 switches off reduced order for solution representation
716 """
717 new_fs=escript.Solution(self.getDomain())
718 if self.getFunctionSpaceForSolution()!=new_fs:
719 self.trace("Full order is used to interpolate solution.")
720 self.__column_function_space=new_fs
721 self.__resetSystem()
722
723 def setReducedOrderForSolutionTo(self,flag=False):
724 """
725 sets order for test functions according to flag
726
727 @param flag: if flag is True, the order reduction is switched on for solution representation, otherwise or
728 if flag is not present order reduction is switched off
729 @type flag: C{bool}
730 """
731 if flag:
732 self.setReducedOrderForSolutionOn()
733 else:
734 self.setReducedOrderForSolutionOff()
735
736 def setReducedOrderForEquationOn(self):
737 """
738 switches on reduced order for equation representation
739 """
740 new_fs=escript.ReducedSolution(self.getDomain())
741 if self.getFunctionSpaceForEquation()!=new_fs:
742 self.trace("Reduced order is used for test functions.")
743 self.__row_function_space=new_fs
744 self.__resetSystem()
745
746 def setReducedOrderForEquationOff(self):
747 """
748 switches off reduced order for equation representation
749 """
750 new_fs=escript.Solution(self.getDomain())
751 if self.getFunctionSpaceForEquation()!=new_fs:
752 self.trace("Full order is used for test functions.")
753 self.__row_function_space=new_fs
754 self.__resetSystem()
755
756 def setReducedOrderForEquationTo(self,flag=False):
757 """
758 sets order for test functions according to flag
759
760 @param flag: if flag is True, the order reduction is switched on for equation representation, otherwise or
761 if flag is not present order reduction is switched off
762 @type flag: C{bool}
763 """
764 if flag:
765 self.setReducedOrderForEquationOn()
766 else:
767 self.setReducedOrderForEquationOff()
768
769 # =============================================================================
770 # private method:
771 # =============================================================================
772 def __checkMatrixType(self):
773 """
774 reassess the matrix type and, if a new matrix is needed, resets the system.
775 """
776 new_matrix_type=self.getDomain().getSystemMatrixTypeId(self.getSolverMethod(),self.isSymmetric())
777 if not new_matrix_type==self.__matrix_type:
778 self.trace("Matrix type is now %d."%new_matrix_type)
779 self.__matrix_type=new_matrix_type
780 self.__resetSystem()
781 #
782 # rebuild switches :
783 #
784 def __invalidateSolution(self):
785 """
786 indicates the PDE has to be resolved if the solution is requested
787 """
788 if self.__solution_isValid: self.trace("PDE has to be resolved.")
789 self.__solution_isValid=False
790
791 def __invalidateOperator(self):
792 """
793 indicates the operator has to be rebuilt next time it is used
794 """
795 if self.__operator_isValid: self.trace("Operator has to be rebuilt.")
796 self.__invalidateSolution()
797 self.__operator_isValid=False
798
799 def __invalidateRightHandSide(self):
800 """
801 indicates the right hand side has to be rebuild next time it is used
802 """
803 if self.__righthandside_isValid: self.trace("Right hand side has to be rebuilt.")
804 self.__invalidateSolution()
805 self.__righthandside_isValid=False
806
807 def __invalidateSystem(self):
808 """
809 annonced that everthing has to be rebuild:
810 """
811 if self.__righthandside_isValid: self.trace("System has to be rebuilt.")
812 self.__invalidateSolution()
813 self.__invalidateOperator()
814 self.__invalidateRightHandSide()
815
816 def __resetSystem(self):
817 """
818 annonced that everthing has to be rebuild:
819 """
820 self.trace("New System is built from scratch.")
821 self.__operator=escript.Operator()
822 self.__operator_isValid=False
823 self.__righthandside=escript.Data()
824 self.__righthandside_isValid=False
825 self.__solution=escript.Data()
826 self.__solution_isValid=False
827 #
828 # system initialization:
829 #
830 def __getNewOperator(self):
831 """
832 returns an instance of a new operator
833 """
834 self.trace("New operator is allocated.")
835 return self.getDomain().newOperator( \
836 self.getNumEquations(), \
837 self.getFunctionSpaceForEquation(), \
838 self.getNumSolutions(), \
839 self.getFunctionSpaceForSolution(), \
840 self.__matrix_type)
841
842 def __getNewRightHandSide(self):
843 """
844 returns an instance of a new right hand side
845 """
846 self.trace("New right hand side is allocated.")
847 if self.getNumEquations()>1:
848 return escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True)
849 else:
850 return escript.Data(0.,(),self.getFunctionSpaceForEquation(),True)
851
852 def __getNewSolution(self):
853 """
854 returns an instance of a new solution
855 """
856 self.trace("New solution is allocated.")
857 if self.getNumSolutions()>1:
858 return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)
859 else:
860 return escript.Data(0.,(),self.getFunctionSpaceForSolution(),True)
861
862 def __makeFreshSolution(self):
863 """
864 makes sure that the solution is instantiated and returns it initialized by zeros
865 """
866 if self.__solution.isEmpty():
867 self.__solution=self.__getNewSolution()
868 else:
869 self.__solution*=0
870 self.trace("Solution is reset to zero.")
871 return self.__solution
872
873 def __makeFreshRightHandSide(self):
874 """
875 makes sure that the right hand side is instantiated and returns it initialized by zeros
876 """
877 if self.__righthandside.isEmpty():
878 self.__righthandside=self.__getNewRightHandSide()
879 else:
880 self.__righthandside*=0
881 self.trace("Right hand side is reset to zero.")
882 return self.__righthandside
883
884 def __makeFreshOperator(self):
885 """
886 makes sure that the operator is instantiated and returns it initialized by zeros
887 """
888 if self.__operator.isEmpty():
889 self.__operator=self.__getNewOperator()
890 else:
891 self.__operator.setValue(0.)
892 self.trace("Operator reset to zero")
893 return self.__operator
894
895 def __applyConstraint(self):
896 """
897 applies the constraints defined by q and r to the system
898 """
899 if not self.isUsingLumping():
900 q=self.getCoefficientOfGeneralPDE("q")
901 r=self.getCoefficientOfGeneralPDE("r")
902 if not q.isEmpty() and not self.__operator.isEmpty():
903 # q is the row and column mask to indicate where constraints are set:
904 row_q=escript.Data(q,self.getFunctionSpaceForEquation())
905 col_q=escript.Data(q,self.getFunctionSpaceForSolution())
906 u=self.__getNewSolution()
907 if r.isEmpty():
908 r_s=self.__getNewSolution()
909 else:
910 r_s=escript.Data(r,self.getFunctionSpaceForSolution())
911 u.copyWithMask(r_s,col_q)
912 if not self.__righthandside.isEmpty():
913 self.__righthandside-=self.__operator*u
914 self.__righthandside=self.copyConstraint(self.__righthandside)
915 self.__operator.nullifyRowsAndCols(row_q,col_q,1.)
916 # =============================================================================
917 # function giving access to coefficients of the general PDE:
918 # =============================================================================
919 def getCoefficientOfGeneralPDE(self,name):
920 """
921 return the value of the coefficient name of the general PDE.
922
923 @note This method is called by the assembling routine it can be overwritten
924 to map coefficients of a particular PDE to the general PDE.
925
926 @param name: name of the coefficient requested.
927 @type name: C{string}
928 @return : the value of the coefficient name
929 @rtype : L{escript.Data}
930 @raise IllegalCoefficient: if name is not one of coefficients
931 "A", "B", "C", "D", "X", "Y", "d", "y", "d_contact", "y_contact", "r" or "q".
932 """
933 if self.hasCoefficientOfGeneralPDE(name):
934 return self.getCoefficient(name)
935 else:
936 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
937
938 def hasCoefficientOfGeneralPDE(self,name):
939 """
940 checks if name is a the name of a coefficient of the general PDE.
941
942 @param name: name of the coefficient enquired.
943 @type name: C{string}
944 @return : True if name is the name of a coefficient of the general PDE. Otherwise False.
945 @rtype : C{bool}
946
947 """
948 return self.__COEFFICIENTS_OF_GENEARL_PDE.has_key(name)
949
950 def createCoefficientOfGeneralPDE(self,name):
951 """
952 returns a new instance of a coefficient for coefficient name of the general PDE
953
954 @param name: name of the coefficient requested.
955 @type name: C{string}
956 @return : a coefficient name initialized to 0.
957 @rtype : L{escript.Data}
958 @raise IllegalCoefficient: if name is not one of coefficients
959 "A", "B", "C", "D", "X", "Y", "d", "y", "d_contact", "y_contact", "r" or "q".
960 """
961 if self.hasCoefficientOfGeneralPDE(name):
962 return escript.Data(0,self.getShapeOfCoefficientOfGeneralPDE(name),self.getFunctionSpaceForCoefficientOfGeneralPDE(name))
963 else:
964 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
965
966 def getFunctionSpaceForCoefficientOfGeneralPDE(self,name):
967 """
968 return the L{escript.FunctionSpace} to be used for coefficient name of the general PDE
969
970 @param name: name of the coefficient enquired.
971 @type name: C{string}
972 @return : the function space to be used for coefficient name
973 @rtype : L{escript.FunctionSpace}
974 @raise IllegalCoefficient: if name is not one of coefficients
975 "A", "B", "C", "D", "X", "Y", "d", "y", "d_contact", "y_contact", "r" or "q".
976 """
977 if self.hasCoefficientOfGeneralPDE(name):
978 return self.__COEFFICIENTS_OF_GENEARL_PDE[name].getFunctionSpace(self.getDomain())
979 else:
980 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
981
982 def getShapeOfCoefficientOfGeneralPDE(self,name):
983 """
984 return the shape of the coefficient name of the general PDE
985
986 @param name: name of the coefficient enquired.
987 @type name: C{string}
988 @return : the shape of the coefficient name
989 @rtype : C{tuple} of C{int}
990 @raise IllegalCoefficient: if name is not one of coefficients
991 "A", "B", "C", "D", "X", "Y", "d", "y", "d_contact", "y_contact", "r" or "q".
992
993 """
994 if self.hasCoefficientOfGeneralPDE(name):
995 return self.__COEFFICIENTS_OF_GENEARL_PDE[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())
996 else:
997 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
998
999 # =============================================================================
1000 # functions giving access to coefficients of a particular PDE implementation:
1001 # =============================================================================
1002 def getCoefficient(self,name):
1003 """
1004 returns the value of the coefficient name
1005
1006 @param name: name of the coefficient requested.
1007 @type name: C{string}
1008 @return : the value of the coefficient name
1009 @rtype : L{escript.Data}
1010 @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1011 """
1012 if self.hasCoefficient(name):
1013 return self.COEFFICIENTS[name].getValue()
1014 else:
1015 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1016
1017 def hasCoefficient(self,name):
1018 """
1019 return True if name is the name of a coefficient
1020
1021 @param name: name of the coefficient enquired.
1022 @type name: C{string}
1023 @return : True if name is the name of a coefficient of the general PDE. Otherwise False.
1024 @rtype : C{bool}
1025 """
1026 return self.COEFFICIENTS.has_key(name)
1027
1028 def createCoefficient(self, name):
1029 """
1030 create a L{escript.Data} object corresponding to coefficient name
1031
1032 @return : a coefficient name initialized to 0.
1033 @rtype : L{escript.Data}
1034 @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1035 """
1036 if self.hasCoefficient(name):
1037 return escript.Data(0.,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name))
1038 else:
1039 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1040
1041 def getFunctionSpaceForCoefficient(self,name):
1042 """
1043 return the L{escript.FunctionSpace} to be used for coefficient name
1044
1045 @param name: name of the coefficient enquired.
1046 @type name: C{string}
1047 @return : the function space to be used for coefficient name
1048 @rtype : L{escript.FunctionSpace}
1049 @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1050 """
1051 if self.hasCoefficient(name):
1052 return self.COEFFICIENTS[name].getFunctionSpace(self.getDomain())
1053 else:
1054 raise ValueError,"unknown coefficient %s requested"%name
1055
1056 def getShapeOfCoefficient(self,name):
1057 """
1058 return the shape of the coefficient name
1059
1060 @param name: name of the coefficient enquired.
1061 @type name: C{string}
1062 @return : the shape of the coefficient name
1063 @rtype : C{tuple} of C{int}
1064 @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1065 """
1066 if self.hasCoefficient(name):
1067 return self.COEFFICIENTS[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())
1068 else:
1069 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1070
1071 def resetCoefficients(self):
1072 """
1073 resets all coefficients to there default values.
1074 """
1075 for i in self.COEFFICIENTS.iterkeys():
1076 self.COEFFICIENTS[i].resetValue()
1077
1078 def alteredCoefficient(self,name):
1079 """
1080 announce that coefficient name has been changed
1081
1082 @param name: name of the coefficient enquired.
1083 @type name: C{string}
1084 @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1085 """
1086 if self.hasCoefficient(name):
1087 self.trace("Coefficient %s has been altered."%name)
1088 if self.COEFFICIENTS[name].isAlteringOperator(): self.__invalidateOperator()
1089 if self.COEFFICIENTS[name].isAlteringRightHandSide(): self.__invalidateRightHandSide()
1090 else:
1091 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1092
1093 def copyConstraint(self,u):
1094 """
1095 copies the constraint into u and returns u.
1096
1097 @param u: a function of rank 0 is a single PDE is solved and of shape (numSolution,) for a system of PDEs
1098 @type u: L{escript.Data}
1099 @return : the input u modified by the constraints.
1100 @rtype : L{escript.Data}
1101 @warning: u is altered if it has the appropriate L{escript.FunctionSpace}
1102
1103 """
1104 q=self.getCoefficientOfGeneralPDE("q")
1105 r=self.getCoefficientOfGeneralPDE("r")
1106 if not q.isEmpty():
1107 if u.isEmpty(): u=escript.Data(0.,q.getShape(),q.getFunctionSpace())
1108 if r.isEmpty():
1109 r=escript.Data(0,u.getShape(),u.getFunctionSpace())
1110 else:
1111 r=escript.Data(r,u.getFunctionSpace())
1112 u.copyWithMask(r,escript.Data(q,u.getFunctionSpace()))
1113 return u
1114
1115 def setValue(self,**coefficients):
1116 """
1117 sets new values to coefficients
1118
1119 @note This method is called by the assembling routine it can be overwritten
1120 to map coefficients of a particular PDE to the general PDE.
1121
1122 @param name: name of the coefficient requested.
1123 @type name: C{string}
1124 @keyword A: value for coefficient A.
1125 @type A: any type that can be interpreted as L{escript.Data} object on L{escript.Function}.
1126 @keyword B: value for coefficient B
1127 @type B: any type that can be interpreted as L{escript.Data} object on L{escript.Function}.
1128 @keyword C: value for coefficient C
1129 @type C: any type that can be interpreted as L{escript.Data} object on L{escript.Function}.
1130 @keyword D: value for coefficient D
1131 @type D: any type that can be interpreted as L{escript.Data} object on L{escript.Function}.
1132 @keyword X: value for coefficient X
1133 @type X: any type that can be interpreted as L{escript.Data} object on L{escript.Function}.
1134 @keyword Y: value for coefficient Y
1135 @type Y: any type that can be interpreted as L{escript.Data} object on L{escript.Function}.
1136 @keyword d: value for coefficient d
1137 @type d: any type that can be interpreted as L{escript.Data} object on L{escript.FunctionOnBoundary}.
1138 @keyword y: value for coefficient y
1139 @type y: any type that can be interpreted as L{escript.Data} object on L{escript.FunctionOnBoundary}.
1140 @keyword d_contact: value for coefficient d_contact
1141 @type d_contact: any type that can be interpreted as L{escript.Data} object on L{escript.FunctionOnContactOne}.
1142 or L{escript.FunctionOnContactZero}.
1143 @keyword y_contact: value for coefficient y_contact
1144 @type y_contact: any type that can be interpreted as L{escript.Data} object on L{escript.FunctionOnContactOne}.
1145 or L{escript.FunctionOnContactZero}.
1146 @keyword r: values prescribed to the solution at the locations of constraints
1147 @type r: any type that can be interpreted as L{escript.Data} object on L{escript.Solution} or L{escript.ReducedSolution}
1148 depending of reduced order is used for the solution.
1149 @keyword q: mask for location of constraints
1150 @type q: any type that can be interpreted as L{escript.Data} object on L{escript.Solution} or L{escript.ReducedSolution}
1151 depending of reduced order is used for the representation of the equation.
1152 @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1153
1154 """
1155 # check if the coefficients are legal:
1156 for i in coefficients.iterkeys():
1157 if not self.hasCoefficient(i):
1158 raise IllegalCoefficient,"Attempt to set unknown coefficient %s"%i
1159 # if the number of unknowns or equations is still unknown we try to estimate them:
1160 if self.__numEquations==None or self.__numSolutions==None:
1161 for i,d in coefficients.iteritems():
1162 if hasattr(d,"shape"):
1163 s=d.shape
1164 elif hasattr(d,"getShape"):
1165 s=d.getShape()
1166 else:
1167 s=numarray.array(d).shape
1168 if s!=None:
1169 # get number of equations and number of unknowns:
1170 res=self.COEFFICIENTS[i].estimateNumEquationsAndNumSolutions(self.getDomain(),s)
1171 if res==None:
1172 raise IllegalCoefficientValue,"Illegal shape %s of coefficient %s"%(s,i)
1173 else:
1174 if self.__numEquations==None: self.__numEquations=res[0]
1175 if self.__numSolutions==None: self.__numSolutions=res[1]
1176 if self.__numEquations==None: raise UndefinedPDEError,"unidententified number of equations"
1177 if self.__numSolutions==None: raise UndefinedPDEError,"unidententified number of solutions"
1178 # now we check the shape of the coefficient if numEquations and numSolutions are set:
1179 for i,d in coefficients.iteritems():
1180 try:
1181 self.COEFFICIENTS[i].setValue(self.getDomain(),self.getNumEquations(),self.getNumSolutions(),d)
1182 except IllegalCoefficientValue,m:
1183 raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
1184 self.alteredCoefficient(i)
1185
1186 # check if the systrem is inhomogeneous:
1187 if len(coefficients)>0 and not self.isUsingLumping():
1188 q=self.getCoefficientOfGeneralPDE("q")
1189 r=self.getCoefficientOfGeneralPDE("r")
1190 homogeneous_constraint=True
1191 if not q.isEmpty() and not r.isEmpty():
1192 if util.Lsup(q*r)>=1.e-13*util.Lsup(r):
1193 self.trace("Inhomogeneous constraint detected.")
1194 self.__invalidateSystem()
1195
1196
1197 def getSystem(self):
1198 """
1199 return the operator and right hand side of the PDE
1200 """
1201 if not self.__operator_isValid or not self.__righthandside_isValid:
1202 if self.isUsingLumping():
1203 if not self.__operator_isValid:
1204 if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution(): raise TypeError,"Lumped matrix requires same order for equations and unknowns"
1205 if not self.getCoefficientOfGeneralPDE("A").isEmpty(): raise Warning,"Lumped matrix does not allow coefficient A"
1206 if not self.getCoefficientOfGeneralPDE("B").isEmpty(): raise Warning,"Lumped matrix does not allow coefficient B"
1207 if not self.getCoefficientOfGeneralPDE("C").isEmpty(): raise Warning,"Lumped matrix does not allow coefficient C"
1208 mat=self.__getNewOperator()
1209 self.getDomain().addPDEToSystem(mat,escript.Data(), \
1210 self.getCoefficientOfGeneralPDE("A"), \
1211 self.getCoefficientOfGeneralPDE("B"), \
1212 self.getCoefficientOfGeneralPDE("C"), \
1213 self.getCoefficientOfGeneralPDE("D"), \
1214 escript.Data(), \
1215 escript.Data(), \
1216 self.getCoefficientOfGeneralPDE("d"), \
1217 escript.Data(),\
1218 self.getCoefficientOfGeneralPDE("d_contact"), \
1219 escript.Data())
1220 self.__operator=1./(mat*escript.Data(1,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True))
1221 del mat
1222 self.trace("New lumped operator has been built.")
1223 self.__operator_isValid=True
1224 if not self.__righthandside_isValid:
1225 self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \
1226 self.getCoefficientOfGeneralPDE("X"), \
1227 self.getCoefficientOfGeneralPDE("Y"),\
1228 self.getCoefficientOfGeneralPDE("y"),\
1229 self.getCoefficientOfGeneralPDE("y_contact"))
1230 self.trace("New right hand side as been built.")
1231 self.__righthandside_isValid=True
1232 else:
1233 if not self.__operator_isValid and not self.__righthandside_isValid:
1234 self.getDomain().addPDEToSystem(self.__makeFreshOperator(),self.__makeFreshRightHandSide(), \
1235 self.getCoefficientOfGeneralPDE("A"), \
1236 self.getCoefficientOfGeneralPDE("B"), \
1237 self.getCoefficientOfGeneralPDE("C"), \
1238 self.getCoefficientOfGeneralPDE("D"), \
1239 self.getCoefficientOfGeneralPDE("X"), \
1240 self.getCoefficientOfGeneralPDE("Y"), \
1241 self.getCoefficientOfGeneralPDE("d"), \
1242 self.getCoefficientOfGeneralPDE("y"), \
1243 self.getCoefficientOfGeneralPDE("d_contact"), \
1244 self.getCoefficientOfGeneralPDE("y_contact"))
1245 self.__applyConstraint()
1246 self.__righthandside=self.copyConstraint(self.__righthandside)
1247 self.trace("New system has been built.")
1248 self.__operator_isValid=True
1249 self.__righthandside_isValid=True
1250 elif not self.__righthandside_isValid:
1251 self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \
1252 self.getCoefficientOfGeneralPDE("X"), \
1253 self.getCoefficientOfGeneralPDE("Y"),\
1254 self.getCoefficientOfGeneralPDE("y"),\
1255 self.getCoefficientOfGeneralPDE("y_contact"))
1256 self.__righthandside=self.copyConstraint(self.__righthandside)
1257 self.trace("New right hand side has been built.")
1258 self.__righthandside_isValid=True
1259 elif not self.__operator_isValid:
1260 self.getDomain().addPDEToSystem(self.__makeFreshOperator(),escript.Data(), \
1261 self.getCoefficientOfGeneralPDE("A"), \
1262 self.getCoefficientOfGeneralPDE("B"), \
1263 self.getCoefficientOfGeneralPDE("C"), \
1264 self.getCoefficientOfGeneralPDE("D"), \
1265 escript.Data(), \
1266 escript.Data(), \
1267 self.getCoefficientOfGeneralPDE("d"), \
1268 escript.Data(),\
1269 self.getCoefficientOfGeneralPDE("d_contact"), \
1270 escript.Data())
1271 self.__applyConstraint()
1272 self.trace("New operator has been built.")
1273 self.__operator_isValid=True
1274 return (self.__operator,self.__righthandside)
1275
1276
1277
1278
1279 class AdvectivePDE(LinearPDE):
1280 """
1281 Class to handle a linear PDE dominated by advective terms:
1282
1283 class to define a linear PDE of the form
1284
1285 \f[
1286 -(A_{ijkl}u_{k,l})_{,j} -(B_{ijk}u_k)_{,j} + C_{ikl}u_{k,l} +D_{ik}u_k = - (X_{ij})_{,j} + Y_i
1287 \f]
1288
1289 with boundary conditons:
1290
1291 \f[
1292 n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_{ik}u_k = - n_j*X_{ij} + y_i
1293 \f]
1294
1295 and contact conditions
1296
1297 \f[
1298 n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d^{contact}_{ik}[u_k] = - n_j*X_{ij} + y^{contact}_{i}
1299 \f]
1300
1301 and constraints:
1302
1303 \f[
1304 u_i=r_i \quad \mathrm{where} \quad q_i>0
1305 \f]
1306 """
1307 def __init__(self,domain,numEquations=0,numSolutions=0,xi=None,debug=False):
1308 LinearPDE.__init__(self,domain,numEquations,numSolutions,debug)
1309 if xi==None:
1310 self.__xi=AdvectivePDE.ELMAN_RAMAGE
1311 else:
1312 self.__xi=xi
1313 self.__Xi=escript.Data()
1314
1315 def __calculateXi(self,peclet_factor,Z,h):
1316 Z_max=util.Lsup(Z)
1317 if Z_max>0.:
1318 return h*self.__xi(Z*peclet_factor)/(Z+Z_max*self.TOL)
1319 else:
1320 return 0.
1321
1322 def setValue(self,**args):
1323 if "A" in args.keys() or "B" in args.keys() or "C" in args.keys(): self.__Xi=escript.Data()
1324 LinearPDE.setValue(**args)
1325
1326 def ELMAN_RAMAGE(P):
1327 """ """
1328 return (P-1.).wherePositive()*0.5*(1.-1./(P+1.e-15))
1329 def SIMPLIFIED_BROOK_HUGHES(P):
1330 """ """
1331 c=(P-3.).whereNegative()
1332 return P/6.*c+1./2.*(1.-c)
1333
1334 def HALF(P):
1335 """ """
1336 return escript.Scalar(0.5,P.getFunctionSpace())
1337
1338 def getXi(self):
1339 if self.__Xi.isEmpty():
1340 B=self.getCoefficient("B")
1341 C=self.getCoefficient("C")
1342 A=self.getCoefficient("A")
1343 h=self.getDomain().getSize()
1344 self.__Xi=escript.Scalar(0.,self.getFunctionSpaceForCoefficient("A"))
1345 if not C.isEmpty() or not B.isEmpty():
1346 if not C.isEmpty() and not B.isEmpty():
1347 Z2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))
1348 if self.getNumEquations()>1:
1349 if self.getNumSolutions()>1:
1350 for i in range(self.getNumEquations()):
1351 for k in range(self.getNumSolutions()):
1352 for l in range(self.getDim()): Z2+=(C[i,k,l]-B[i,l,k])**2
1353 else:
1354 for i in range(self.getNumEquations()):
1355 for l in range(self.getDim()): Z2+=(C[i,l]-B[i,l])**2
1356 else:
1357 if self.getNumSolutions()>1:
1358 for k in range(self.getNumSolutions()):
1359 for l in range(self.getDim()): Z2+=(C[k,l]-B[l,k])**2
1360 else:
1361 for l in range(self.getDim()): Z2+=(C[l]-B[l])**2
1362 length_of_Z=util.sqrt(Z2)
1363 elif C.isEmpty():
1364 length_of_Z=util.length(B)
1365 else:
1366 length_of_Z=util.length(C)
1367
1368 Z_max=util.Lsup(length_of_Z)
1369 if Z_max>0.:
1370 length_of_A=util.length(A)
1371 A_max=util.Lsup(length_of_A)
1372 if A_max>0:
1373 inv_A=1./(length_of_A+A_max*self.TOL)
1374 else:
1375 inv_A=1./self.TOL
1376 peclet_number=length_of_Z*h/2*inv_A
1377 xi=self.__xi(peclet_number)
1378 self.__Xi=h*xi/(length_of_Z+Z_max*self.TOL)
1379 print "@ preclet number = %e"%util.Lsup(peclet_number),util.Lsup(xi),util.Lsup(length_of_Z)
1380 return self.__Xi
1381
1382
1383 def getCoefficientOfGeneralPDE(self,name):
1384 """
1385 return the value of the coefficient name of the general PDE
1386
1387 @param name:
1388 """
1389 if not self.getNumEquations() == self.getNumSolutions():
1390 raise ValueError,"AdvectivePDE expects the number of solution componets and the number of equations to be equal."
1391
1392 if name == "A" :
1393 A=self.getCoefficient("A")
1394 B=self.getCoefficient("B")
1395 C=self.getCoefficient("C")
1396 if B.isEmpty() and C.isEmpty():
1397 Aout=A
1398 else:
1399 if A.isEmpty():
1400 Aout=self.createNewCoefficient("A")
1401 else:
1402 Aout=A[:]
1403 Xi=self.getXi()
1404 if self.getNumEquations()>1:
1405 for i in range(self.getNumEquations()):
1406 for j in range(self.getDim()):
1407 for k in range(self.getNumSolutions()):
1408 for l in range(self.getDim()):
1409 if not C.isEmpty() and not B.isEmpty():
1410 for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*(C[p,i,j]-B[p,j,i])*(C[p,k,l]-B[p,l,k])
1411 elif C.isEmpty():
1412 for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*B[p,j,i]*B[p,l,k]
1413 else:
1414 for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*C[p,i,j]*C[p,k,l]
1415 else:
1416 for j in range(self.getDim()):
1417 for l in range(self.getDim()):
1418 if not C.isEmpty() and not B.isEmpty():
1419 Aout[j,l]+=Xi*(C[j]-B[j])*(C[l]-B[l])
1420 elif C.isEmpty():
1421 Aout[j,l]+=Xi*B[j]*B[l]
1422 else:
1423 Aout[j,l]+=Xi*C[j]*C[l]
1424 return Aout
1425 elif name == "B" :
1426 B=self.getCoefficient("B")
1427 C=self.getCoefficient("C")
1428 D=self.getCoefficient("D")
1429 if C.isEmpty() or D.isEmpty():
1430 Bout=B
1431 else:
1432 Xi=self.getXi()
1433 if B.isEmpty():
1434 Bout=self.createNewCoefficient("B")
1435 else:
1436 Bout=B[:]
1437 if self.getNumEquations()>1:
1438 for k in range(self.getNumSolutions()):
1439 for p in range(self.getNumEquations()):
1440 tmp=Xi*D[p,k]
1441 for i in range(self.getNumEquations()):
1442 for j in range(self.getDim()):
1443 Bout[i,j,k]+=tmp*C[p,i,j]
1444 else:
1445 tmp=Xi*D
1446 for j in range(self.getDim()): Bout[j]+=tmp*C[j]
1447 return Bout
1448 elif name == "C" :
1449 B=self.getCoefficient("B")
1450 C=self.getCoefficient("C")
1451 D=self.getCoefficient("D")
1452 if B.isEmpty() or D.isEmpty():
1453 Cout=C
1454 else:
1455 Xi=self.getXi()
1456 if C.isEmpty():
1457 Cout=self.createNewCoefficient("C")
1458 else:
1459 Cout=C[:]
1460 if self.getNumEquations()>1:
1461 for k in range(self.getNumSolutions()):
1462 for p in range(self.getNumEquations()):
1463 tmp=Xi*D[p,k]
1464 for i in range(self.getNumEquations()):
1465 for l in range(self.getDim()):
1466 Cout[i,k,l]+=tmp*B[p,l,i]
1467 else:
1468 tmp=Xi*D
1469 for j in range(self.getDim()): Cout[j]+=tmp*B[j]
1470 return Cout
1471 elif name == "D" :
1472 return self.getCoefficient("D")
1473 elif name == "X" :
1474 X=self.getCoefficient("X")
1475 Y=self.getCoefficient("Y")
1476 B=self.getCoefficient("B")
1477 C=self.getCoefficient("C")
1478 if Y.isEmpty() or (B.isEmpty() and C.isEmpty()):
1479 Xout=X
1480 else:
1481 if X.isEmpty():
1482 Xout=self.createNewCoefficient("X")
1483 else:
1484 Xout=X[:]
1485 Xi=self.getXi()
1486 if self.getNumEquations()>1:
1487 for p in range(self.getNumEquations()):
1488 tmp=Xi*Y[p]
1489 for i in range(self.getNumEquations()):
1490 for j in range(self.getDim()):
1491 if not C.isEmpty() and not B.isEmpty():
1492 Xout[i,j]+=tmp*(C[p,i,j]-B[p,j,i])
1493 elif C.isEmpty():
1494 Xout[i,j]-=tmp*B[p,j,i]
1495 else:
1496 Xout[i,j]+=tmp*C[p,i,j]
1497 else:
1498 tmp=Xi*Y
1499 for j in range(self.getDim()):
1500 if not C.isEmpty() and not B.isEmpty():
1501 Xout[j]+=tmp*(C[j]-B[j])
1502 elif C.isEmpty():
1503 Xout[j]-=tmp*B[j]
1504 else:
1505 Xout[j]+=tmp*C[j]
1506 return Xout
1507 elif name == "Y" :
1508 return self.getCoefficient("Y")
1509 elif name == "d" :
1510 return self.getCoefficient("d")
1511 elif name == "y" :
1512 return self.getCoefficient("y")
1513 elif name == "d_contact" :
1514 return self.getCoefficient("d_contact")
1515 elif name == "y_contact" :
1516 return self.getCoefficient("y_contact")
1517 elif name == "r" :
1518 return self.getCoefficient("r")
1519 elif name == "q" :
1520 return self.getCoefficient("q")
1521 else:
1522 raise SystemError,"unknown PDE coefficient %s",name
1523
1524
1525 class Poisson(LinearPDE):
1526 """
1527 Class to define a Poisson equation problem:
1528
1529 class to define a linear PDE of the form
1530 \f[
1531 -u_{,jj} = f
1532 \f]
1533
1534 with boundary conditons:
1535
1536 \f[
1537 n_j*u_{,j} = 0
1538 \f]
1539
1540 and constraints:
1541
1542 \f[
1543 u=0 \quad \mathrm{where} \quad q>0
1544 \f]
1545 """
1546
1547 def __init__(self,domain,f=escript.Data(),q=escript.Data(),debug=False):
1548 LinearPDE.__init__(self,domain,1,1,debug)
1549 self.COEFFICIENTS={"f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1550 "q": PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.BOTH)}
1551 self.setSymmetryOn()
1552 self.setValue(f,q)
1553
1554 def setValue(self,f=escript.Data(),q=escript.Data()):
1555 """set value of PDE parameters f and q"""
1556 self._LinearPDE__setValue(f=f,q=q)
1557
1558 def getCoefficientOfGeneralPDE(self,name):
1559 """
1560 return the value of the coefficient name of the general PDE
1561
1562 @param name:
1563 """
1564 if name == "A" :
1565 return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))
1566 elif name == "B" :
1567 return escript.Data()
1568 elif name == "C" :
1569 return escript.Data()
1570 elif name == "D" :
1571 return escript.Data()
1572 elif name == "X" :
1573 return escript.Data()
1574 elif name == "Y" :
1575 return self.getCoefficient("f")
1576 elif name == "d" :
1577 return escript.Data()
1578 elif name == "y" :
1579 return escript.Data()
1580 elif name == "d_contact" :
1581 return escript.Data()
1582 elif name == "y_contact" :
1583 return escript.Data()
1584 elif name == "r" :
1585 return escript.Data()
1586 elif name == "q" :
1587 return self.getCoefficient("q")
1588 else:
1589 raise SystemError,"unknown PDE coefficient %s",name
1590
1591 class LameEquation(LinearPDE):
1592 """
1593 Class to define a Lame equation problem:
1594
1595 class to define a linear PDE of the form
1596 \f[
1597 -(\mu (u_{i,j}+u_{j,i}))_{,j} - \lambda u_{j,ji}} = F_i -\sigma_{ij,j}
1598 \f]
1599
1600 with boundary conditons:
1601
1602 \f[
1603 n_j(\mu (u_{i,j}+u_{j,i})-sigma_{ij}) + n_i\lambda u_{j,j} = f_i
1604 \f]
1605
1606 and constraints:
1607
1608 \f[
1609 u_i=r_i \quad \mathrm{where} \quad q_i>0
1610 \f]
1611 """
1612
1613 def __init__(self,domain,f=escript.Data(),q=escript.Data(),debug=False):
1614 LinearPDE.__init__(self,domain,domain.getDim(),domain.getDim(),debug)
1615 self.COEFFICIENTS={ "lame_lambda" : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),
1616 "lame_mu" : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),
1617 "F" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1618 "sigma" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM),PDECoefficient.RIGHTHANDSIDE),
1619 "f" : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1620 "r" : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.BOTH),
1621 "q" : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.BOTH)}
1622 self.setSymmetryOn()
1623
1624 def setValue(self,lame_lambda=escript.Data(),lame_mu=escript.Data(),F=escript.Data(),sigma=escript.Data(),f=escript.Data(),r=escript.Data(),q=escript.Data()):
1625 """set value of PDE parameters"""
1626 self._LinearPDE__setValue(lame_lambda=lame_lambda, \
1627 lame_mu=lame_mu, \
1628 F=F, \
1629 sigma=sigma, \
1630 f=f, \
1631 r=r, \
1632 q=q)
1633 def getCoefficientOfGeneralPDE(self,name):
1634 """
1635 return the value of the coefficient name of the general PDE
1636
1637 @param name:
1638 """
1639 if name == "A" :
1640 out =self.createCoefficientOfGeneralPDE("A")
1641 for i in range(self.getDim()):
1642 for j in range(self.getDim()):
1643 out[i,i,j,j] += self.getCoefficient("lame_lambda")
1644 out[i,j,j,i] += self.getCoefficient("lame_mu")
1645 out[i,j,i,j] += self.getCoefficient("lame_mu")
1646 return out
1647 elif name == "B" :
1648 return escript.Data()
1649 elif name == "C" :
1650 return escript.Data()
1651 elif name == "D" :
1652 return escript.Data()
1653 elif name == "X" :
1654 return self.getCoefficient("sigma")
1655 elif name == "Y" :
1656 return self.getCoefficient("F")
1657 elif name == "d" :
1658 return escript.Data()
1659 elif name == "y" :
1660 return self.getCoefficient("f")
1661 elif name == "d_contact" :
1662 return escript.Data()
1663 elif name == "y_contact" :
1664 return escript.Data()
1665 elif name == "r" :
1666 return self.getCoefficient("r")
1667 elif name == "q" :
1668 return self.getCoefficient("q")
1669 else:
1670 raise SystemError,"unknown PDE coefficient %s",name
1671
1672 # $Log$
1673 # Revision 1.11 2005/08/23 01:24:28 jgs
1674 # Merge of development branch dev-02 back to main trunk on 2005-08-23
1675 #
1676 # Revision 1.10 2005/08/12 01:45:36 jgs
1677 # erge of development branch dev-02 back to main trunk on 2005-08-12
1678 #
1679 # Revision 1.9.2.4 2005/08/22 07:11:09 gross
1680 # some problems with LinearPDEs fixed.
1681 #
1682 # Revision 1.9.2.3 2005/08/18 04:48:48 gross
1683 # the methods SetLumping*() are removed. Lumping is set trough setSolverMethod(LinearPDE.LUMPING)
1684 #
1685 # Revision 1.9.2.2 2005/08/18 04:39:32 gross
1686 # the constants have been removed from util.py as they not needed anymore. PDE related constants are accessed through LinearPDE attributes now
1687 #
1688 # Revision 1.9.2.1 2005/07/29 07:10:27 gross
1689 # new functions in util and a new pde type in linearPDEs
1690 #
1691 # Revision 1.1.2.25 2005/07/28 04:21:09 gross
1692 # Lame equation: (linear elastic, isotropic) added
1693 #
1694 # Revision 1.1.2.24 2005/07/22 06:37:11 gross
1695 # some extensions to modellib and linearPDEs
1696 #
1697 # Revision 1.1.2.23 2005/05/13 00:55:20 cochrane
1698 # Fixed up some docstrings. Moved module-level functions to top of file so
1699 # that epydoc and doxygen can pick them up properly.
1700 #
1701 # Revision 1.1.2.22 2005/05/12 11:41:30 gross
1702 # some basic Models have been added
1703 #
1704 # Revision 1.1.2.21 2005/05/12 07:16:12 cochrane
1705 # Moved ELMAN_RAMAGE, SIMPLIFIED_BROOK_HUGHES, and HALF functions to bottom of
1706 # file so that the AdvectivePDE class is picked up by doxygen. Some
1707 # reformatting of docstrings. Addition of code to make equations come out
1708 # as proper LaTeX.
1709 #
1710 # Revision 1.1.2.20 2005/04/15 07:09:08 gross
1711 # some problems with functionspace and linearPDEs fixed.
1712 #
1713 # Revision 1.1.2.19 2005/03/04 05:27:07 gross
1714 # bug in SystemPattern fixed.
1715 #
1716 # Revision 1.1.2.18 2005/02/08 06:16:45 gross
1717 # Bugs in AdvectivePDE fixed, AdvectiveTest is stable but more testing is needed
1718 #
1719 # Revision 1.1.2.17 2005/02/08 05:56:19 gross
1720 # Reference Number handling added
1721 #
1722 # Revision 1.1.2.16 2005/02/07 04:41:28 gross
1723 # some function exposed to python to make mesh merging running
1724 #
1725 # Revision 1.1.2.15 2005/02/03 00:14:44 gross
1726 # timeseries add and ESySParameter.py renames esysXML.py for consistence
1727 #
1728 # Revision 1.1.2.14 2005/02/01 06:44:10 gross
1729 # new implementation of AdvectivePDE which now also updates right hand side. systems of PDEs are still not working
1730 #
1731 # Revision 1.1.2.13 2005/01/25 00:47:07 gross
1732 # updates in the documentation
1733 #
1734 # Revision 1.1.2.12 2005/01/12 01:28:04 matt
1735 # Added createCoefficient method for linearPDEs.
1736 #
1737 # Revision 1.1.2.11 2005/01/11 01:55:34 gross
1738 # a problem in linearPDE class fixed
1739 #
1740 # Revision 1.1.2.10 2005/01/07 01:13:29 gross
1741 # some bugs in linearPDE fixed
1742 #
1743 # Revision 1.1.2.9 2005/01/06 06:24:58 gross
1744 # some bugs in slicing fixed
1745 #
1746 # Revision 1.1.2.8 2005/01/05 04:21:40 gross
1747 # FunctionSpace checking/matchig in slicing added
1748 #
1749 # Revision 1.1.2.7 2004/12/29 10:03:41 gross
1750 # bug in setValue fixed
1751 #
1752 # Revision 1.1.2.6 2004/12/29 05:29:59 gross
1753 # AdvectivePDE successfully tested for Peclet number 1000000. there is still a problem with setValue and Data()
1754 #
1755 # Revision 1.1.2.5 2004/12/29 00:18:41 gross
1756 # AdvectivePDE added
1757 #
1758 # Revision 1.1.2.4 2004/12/24 06:05:41 gross
1759 # some changes in linearPDEs to add AdevectivePDE
1760 #
1761 # Revision 1.1.2.3 2004/12/16 00:12:34 gross
1762 # __init__ of LinearPDE does not accept any coefficient anymore
1763 #
1764 # Revision 1.1.2.2 2004/12/14 03:55:01 jgs
1765 # *** empty log message ***
1766 #
1767 # Revision 1.1.2.1 2004/12/12 22:53:47 gross
1768 # linearPDE has been renamed LinearPDE
1769 #
1770 # Revision 1.1.1.1.2.7 2004/12/07 10:13:08 gross
1771 # GMRES added
1772 #
1773 # Revision 1.1.1.1.2.6 2004/12/07 03:19:50 gross
1774 # options for GMRES and PRES20 added
1775 #
1776 # Revision 1.1.1.1.2.5 2004/12/01 06:25:15 gross
1777 # some small changes
1778 #
1779 # Revision 1.1.1.1.2.4 2004/11/24 01:50:21 gross
1780 # Finley solves 4M unknowns now
1781 #
1782 # Revision 1.1.1.1.2.3 2004/11/15 06:05:26 gross
1783 # poisson solver added
1784 #
1785 # Revision 1.1.1.1.2.2 2004/11/12 06:58:15 gross
1786 # a lot of changes to get the linearPDE class running: most important change is that there is no matrix format exposed to the user anymore. the format is chosen by the Domain according to the solver and symmetry
1787 #
1788 # Revision 1.1.1.1.2.1 2004/10/28 22:59:22 gross
1789 # finley's RecTest.py is running now: problem in SystemMatrixAdapater fixed
1790 #
1791 # Revision 1.1.1.1 2004/10/26 06:53:56 jgs
1792 # initial import of project esys2
1793 #
1794 # Revision 1.3.2.3 2004/10/26 06:43:48 jgs
1795 # committing Lutz's and Paul's changes to brach jgs
1796 #
1797 # Revision 1.3.4.1 2004/10/20 05:32:51 cochrane
1798 # Added incomplete Doxygen comments to files, or merely put the docstrings that already exist into Doxygen form.
1799 #
1800 # Revision 1.3 2004/09/23 00:53:23 jgs
1801 # minor fixes
1802 #
1803 # Revision 1.1 2004/08/28 12:58:06 gross
1804 # SimpleSolve is not running yet: problem with == of functionsspace
1805 #
1806 #

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26