/[escript]/trunk/esys2/escript/py_src/linearPDEs.py
ViewVC logotype

Contents of /trunk/esys2/escript/py_src/linearPDEs.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 110 - (show annotations)
Mon Feb 14 04:14:42 2005 UTC (14 years, 2 months ago) by jgs
File MIME type: text/x-python
File size: 47700 byte(s)
*** empty log message ***

1 # $Id$
2
3 ## @file linearPDEs.py
4
5 """
6 @brief Functions and classes for linear PDEs
7 """
8
9 import escript
10 import util
11 import numarray
12
13
14 def _CompTuple2(t1,t2):
15 """
16 @brief
17
18 @param t1
19 @param t2
20 """
21 dif=t1[0]+t1[1]-(t2[0]+t2[1])
22 if dif<0: return 1
23 elif dif>0: return -1
24 else: return 0
25
26 class PDECoefficient:
27 """
28 @brief
29 """
30 # identifier for location of Data objects defining COEFFICIENTS
31 INTERIOR=0
32 BOUNDARY=1
33 CONTACT=2
34 CONTINUOUS=3
35 # identifier in the pattern of COEFFICIENTS:
36 # the pattern is a tuple of EQUATION,SOLUTION,DIM where DIM represents the spatial dimension, EQUATION the number of equations and SOLUTION the
37 # number of unknowns.
38 EQUATION=3
39 SOLUTION=4
40 DIM=5
41 # indicator for what is altered if the coefficient is altered:
42 OPERATOR=5
43 RIGHTHANDSIDE=6
44 BOTH=7
45 def __init__(self,where,pattern,altering):
46 """
47 @brief Initialise a PDE Coefficient type
48 """
49 self.what=where
50 self.pattern=pattern
51 self.altering=altering
52 self.resetValue()
53
54 def resetValue(self):
55 """
56 @brief resets coefficient value to default
57 """
58 self.value=escript.Data()
59
60 def getFunctionSpace(self,domain):
61 """
62 @brief defines the FunctionSpace of the coefficient on the domain
63
64 @param domain
65 """
66 if self.what==self.INTERIOR: return escript.Function(domain)
67 elif self.what==self.BOUNDARY: return escript.FunctionOnBoundary(domain)
68 elif self.what==self.CONTACT: return escript.FunctionOnContactZero(domain)
69 elif self.what==self.CONTINUOUS: return escript.ContinuousFunction(domain)
70
71 def getValue(self):
72 """
73 @brief returns the value of the coefficient:
74 """
75 return self.value
76
77 def setValue(self,newValue):
78 """
79 @brief set the value of the coefficient to new value
80 """
81 self.value=newValue
82
83 def isAlteringOperator(self):
84 """
85 @brief return true if the operator of the PDE is changed when the coefficient is changed
86 """
87 if self.altering==self.OPERATOR or self.altering==self.BOTH:
88 return not None
89 else:
90 return None
91
92 def isAlteringRightHandSide(self):
93 """
94 @brief return true if the right hand side of the PDE is changed when the coefficient is changed
95 """
96 if self.altering==self.RIGHTHANDSIDE or self.altering==self.BOTH:
97 return not None
98 else:
99 return None
100
101 def estimateNumEquationsAndNumSolutions(self,shape=(),dim=3):
102 """
103 @brief tries to estimate the number of equations in a given tensor shape for a given spatial dimension dim
104
105 @param shape
106 @param dim
107 """
108 if len(shape)>0:
109 num=max(shape)+1
110 else:
111 num=1
112 search=[]
113 for u in range(num):
114 for e in range(num):
115 search.append((e,u))
116 search.sort(_CompTuple2)
117 for item in search:
118 s=self.buildShape(item[0],item[1],dim)
119 if len(s)==0 and len(shape)==0:
120 return (1,1)
121 else:
122 if s==shape: return item
123 return None
124
125 def buildShape(self,e=1,u=1,dim=3):
126 """
127 @brief builds the required shape for a given number of equations e, number of unknowns u and spatial dimension dim
128
129 @param e
130 @param u
131 @param dim
132 """
133 s=()
134 for i in self.pattern:
135 if i==self.EQUATION:
136 if e>1: s=s+(e,)
137 elif i==self.SOLUTION:
138 if u>1: s=s+(u,)
139 else:
140 s=s+(dim,)
141 return s
142
143 class LinearPDE:
144 """
145 @brief Class to handel a linear PDE
146
147 class to define a linear PDE of the form
148
149 -(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
150
151 with boundary conditons:
152
153 n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_{ik}u_k = - n_j*X_{ij} + y_i
154
155 and contact conditions
156
157 n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_contact_{ik}[u_k] = - n_j*X_{ij} + y_contact_i
158
159 and constraints:
160
161 u_i=r_i where q_i>0
162
163 """
164 TOL=1.e-13
165 DEFAULT_METHOD=util.DEFAULT_METHOD
166 DIRECT=util.DIRECT
167 CHOLEVSKY=util.CHOLEVSKY
168 PCG=util.PCG
169 CR=util.CR
170 CGS=util.CGS
171 BICGSTAB=util.BICGSTAB
172 SSOR=util.SSOR
173 GMRES=util.GMRES
174 PRES20=util.PRES20
175
176 def __init__(self,domain,numEquations=0,numSolutions=0):
177 """
178 @brief initializes a new linear PDE.
179
180 @param args
181 """
182 # COEFFICIENTS can be overwritten by subclasses:
183 self.COEFFICIENTS={
184 "A" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM,PDECoefficient.SOLUTION,PDECoefficient.DIM),PDECoefficient.OPERATOR),
185 "B" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),
186 "C" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION,PDECoefficient.DIM),PDECoefficient.OPERATOR),
187 "D" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),
188 "X" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM),PDECoefficient.RIGHTHANDSIDE),
189 "Y" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
190 "d" : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),
191 "y" : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
192 "d_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),
193 "y_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
194 "r" : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
195 "q" : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.SOLUTION,),PDECoefficient.BOTH)}
196
197 # initialize attributes
198 self.__debug=None
199 self.__domain=domain
200 self.__numEquations=numEquations
201 self.__numSolutions=numSolutions
202 self.cleanCoefficients()
203
204 self.__operator=escript.Operator()
205 self.__operator_isValid=False
206 self.__righthandside=escript.Data()
207 self.__righthandside_isValid=False
208 self.__solution=escript.Data()
209 self.__solution_isValid=False
210
211 # set some default values:
212 self.__homogeneous_constraint=True
213 self.__row_function_space=escript.Solution(self.__domain)
214 self.__column_function_space=escript.Solution(self.__domain)
215 self.__tolerance=1.e-8
216 self.__solver_method=util.DEFAULT_METHOD
217 self.__matrix_type=self.__domain.getSystemMatrixTypeId(util.DEFAULT_METHOD,False)
218 self.__sym=False
219 self.__lumping=False
220
221 def createCoefficient(self, name):
222 """
223 @brief create a data object corresponding to coefficient name
224 @param name
225 """
226 return escript.Data(shape = getShapeOfCoefficient(name), \
227 what = getFunctionSpaceForCoefficient(name))
228
229 def __del__(self):
230 pass
231
232 def getCoefficient(self,name):
233 """
234 @brief return the value of the parameter name
235
236 @param name
237 """
238 return self.COEFFICIENTS[name].getValue()
239
240 def getCoefficientOfPDE(self,name):
241 """
242 @brief return the value of the coefficient name of the general PDE. This method is called by the assembling routine
243 it can be overwritten to map coefficients of a particualr PDE to the general PDE.
244 @param name
245 """
246 return self.getCoefficient(name)
247
248 def hasCoefficient(self,name):
249 """
250 @brief return true if name is the name of a coefficient
251
252 @param name
253 """
254 return self.COEFFICIENTS.has_key(name)
255
256 def getFunctionSpaceForEquation(self):
257 """
258 @brief return true if the test functions should use reduced order
259 """
260 return self.__row_function_space
261
262 def getFunctionSpaceForSolution(self):
263 """
264 @brief return true if the interpolation of the solution should use reduced order
265 """
266 return self.__column_function_space
267
268 def setValue(self,**coefficients):
269 """
270 @brief sets new values to coefficients
271
272 @param coefficients
273 """
274 self._setValue(**coefficients)
275
276
277 def cleanCoefficients(self):
278 """
279 @brief resets all coefficients to default values.
280 """
281 for i in self.COEFFICIENTS.iterkeys():
282 self.COEFFICIENTS[i].resetValue()
283
284 def createNewCoefficient(self,name):
285 """
286 @brief returns a new coefficient appropriate for coefficient name:
287 """
288 return escript.Data(0,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name))
289
290
291 def getShapeOfCoefficient(self,name):
292 """
293 @brief return the shape of the coefficient name
294
295 @param name
296 """
297 if self.hasCoefficient(name):
298 return self.COEFFICIENTS[name].buildShape(self.getNumEquations(),self.getNumSolutions(),self.getDomain().getDim())
299 else:
300 raise ValueError,"Solution coefficient %s requested"%name
301
302 def getFunctionSpaceForCoefficient(self,name):
303 """
304 @brief return the atoms of the coefficient name
305
306 @param name
307 """
308 if self.hasCoefficient(name):
309 return self.COEFFICIENTS[name].getFunctionSpace(self.getDomain())
310 else:
311 raise ValueError,"Solution coefficient %s requested"%name
312
313 def alteredCoefficient(self,name):
314 """
315 @brief annonced that coefficient name has been changed
316
317 @param name
318 """
319 if self.hasCoefficient(name):
320 if self.COEFFICIENTS[name].isAlteringOperator(): self.__rebuildOperator()
321 if self.COEFFICIENTS[name].isAlteringRightHandSide(): self.__rebuildRightHandSide()
322 else:
323 raise ValueError,"unknown coefficient %s requested"%name
324
325 # ===== debug ==============================================================
326 def setDebugOn(self):
327 """
328 @brief
329 """
330 self.__debug=not None
331
332 def setDebugOff(self):
333 """
334 @brief
335 """
336 self.__debug=None
337
338 def debug(self):
339 """
340 @brief returns true if the PDE is in the debug mode
341 """
342 return self.__debug
343
344 #===== Lumping ===========================
345 def setLumpingOn(self):
346 """
347 @brief indicates to use matrix lumping
348 """
349 if not self.isUsingLumping():
350 if self.debug() : print "PDE Debug: lumping is set on"
351 self.__rebuildOperator()
352 self.__lumping=True
353
354 def setLumpingOff(self):
355 """
356 @brief switches off matrix lumping
357 """
358 if self.isUsingLumping():
359 if self.debug() : print "PDE Debug: lumping is set off"
360 self.__rebuildOperator()
361 self.__lumping=False
362
363 def setLumping(self,flag=False):
364 """
365 @brief set the matrix lumping flag to flag
366 """
367 if flag:
368 self.setLumpingOn()
369 else:
370 self.setLumpingOff()
371
372 def isUsingLumping(self):
373 """
374 @brief
375 """
376 return self.__lumping
377
378 #============ method business =========================================================
379 def setSolverMethod(self,solver=util.DEFAULT_METHOD):
380 """
381 @brief sets a new solver
382 """
383 if not solver==self.getSolverMethod():
384 self.__solver_method=solver
385 if self.debug() : print "PDE Debug: New solver is %s"%solver
386 self.__checkMatrixType()
387
388 def getSolverMethod(self):
389 """
390 @brief returns the solver method
391 """
392 return self.__solver_method
393
394 #============ tolerance business =========================================================
395 def setTolerance(self,tol=1.e-8):
396 """
397 @brief resets the tolerance to tol.
398 """
399 if not tol>0:
400 raise ValueException,"Tolerance as to be positive"
401 if tol<self.getTolerance(): self.__rebuildSolution()
402 if self.debug() : print "PDE Debug: New tolerance %e",tol
403 self.__tolerance=tol
404 return
405 def getTolerance(self):
406 """
407 @brief returns the tolerance set for the solution
408 """
409 return self.__tolerance
410
411 #===== symmetry flag ==========================
412 def isSymmetric(self):
413 """
414 @brief returns true is the operator is considered to be symmetric
415 """
416 return self.__sym
417
418 def setSymmetryOn(self):
419 """
420 @brief sets the symmetry flag to true
421 """
422 if not self.isSymmetric():
423 if self.debug() : print "PDE Debug: Operator is set to be symmetric"
424 self.__sym=True
425 self.__checkMatrixType()
426
427 def setSymmetryOff(self):
428 """
429 @brief sets the symmetry flag to false
430 """
431 if self.isSymmetric():
432 if self.debug() : print "PDE Debug: Operator is set to be unsymmetric"
433 self.__sym=False
434 self.__checkMatrixType()
435
436 def setSymmetryTo(self,flag=False):
437 """
438 @brief sets the symmetry flag to flag
439
440 @param flag
441 """
442 if flag:
443 self.setSymmetryOn()
444 else:
445 self.setSymmetryOff()
446
447 #===== order reduction ==========================
448 def setReducedOrderOn(self):
449 """
450 @brief switches to on reduced order
451 """
452 self.setReducedOrderForSolutionOn()
453 self.setReducedOrderForEquationOn()
454
455 def setReducedOrderOff(self):
456 """
457 @brief switches to full order
458 """
459 self.setReducedOrderForSolutionOff()
460 self.setReducedOrderForEquationOff()
461
462 def setReducedOrderTo(self,flag=False):
463 """
464 @brief sets order according to flag
465
466 @param flag
467 """
468 self.setReducedOrderForSolutionTo(flag)
469 self.setReducedOrderForEquationTo(flag)
470
471
472 #===== order reduction solution ==========================
473 def setReducedOrderForSolutionOn(self):
474 """
475 @brief switches to reduced order to interpolate solution
476 """
477 new_fs=escript.ReducedSolution(self.getDomain())
478 if self.getFunctionSpaceForSolution()!=new_fs:
479 if self.debug() : print "PDE Debug: Reduced order is used to interpolate solution."
480 self.__column_function_space=new_fs
481 self.__rebuildSystem(deep=True)
482
483 def setReducedOrderForSolutionOff(self):
484 """
485 @brief switches to full order to interpolate solution
486 """
487 new_fs=escript.Solution(self.getDomain())
488 if self.getFunctionSpaceForSolution()!=new_fs:
489 if self.debug() : print "PDE Debug: Full order is used to interpolate solution."
490 self.__column_function_space=new_fs
491 self.__rebuildSystem(deep=True)
492
493 def setReducedOrderForSolutionTo(self,flag=False):
494 """
495 @brief sets order for test functions according to flag
496
497 @param flag
498 """
499 if flag:
500 self.setReducedOrderForSolutionOn()
501 else:
502 self.setReducedOrderForSolutionOff()
503
504 #===== order reduction equation ==========================
505 def setReducedOrderForEquationOn(self):
506 """
507 @brief switches to reduced order for test functions
508 """
509 new_fs=escript.ReducedSolution(self.getDomain())
510 if self.getFunctionSpaceForEquation()!=new_fs:
511 if self.debug() : print "PDE Debug: Reduced order is used for test functions."
512 self.__row_function_space=new_fs
513 self.__rebuildSystem(deep=True)
514
515 def setReducedOrderForEquationOff(self):
516 """
517 @brief switches to full order for test functions
518 """
519 new_fs=escript.Solution(self.getDomain())
520 if self.getFunctionSpaceForEquation()!=new_fs:
521 if self.debug() : print "PDE Debug: Full order is used for test functions."
522 self.__row_function_space=new_fs
523 self.__rebuildSystem(deep=True)
524
525 def setReducedOrderForEquationTo(self,flag=False):
526 """
527 @brief sets order for test functions according to flag
528
529 @param flag
530 """
531 if flag:
532 self.setReducedOrderForEquationOn()
533 else:
534 self.setReducedOrderForEquationOff()
535
536 # ==== initialization =====================================================================
537 def __makeNewOperator(self):
538 """
539 @brief
540 """
541 return self.getDomain().newOperator( \
542 self.getNumEquations(), \
543 self.getFunctionSpaceForEquation(), \
544 self.getNumSolutions(), \
545 self.getFunctionSpaceForSolution(), \
546 self.__matrix_type)
547
548 def __makeNewRightHandSide(self):
549 """
550 @brief
551 """
552 return escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True)
553
554 def __makeNewSolution(self):
555 """
556 @brief
557 """
558 return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)
559
560 def __getFreshOperator(self):
561 """
562 @brief
563 """
564 if self.__operator.isEmpty():
565 self.__operator=self.__makeNewOperator()
566 if self.debug() : print "PDE Debug: New operator allocated"
567 else:
568 self.__operator.setValue(0.)
569 self.__operator.resetSolver()
570 if self.debug() : print "PDE Debug: Operator reset to zero"
571 return self.__operator
572
573 def __getFreshRightHandSide(self):
574 """
575 @brief
576 """
577 if self.__righthandside.isEmpty():
578 self.__righthandside=self.__makeNewRightHandSide()
579 if self.debug() : print "PDE Debug: New right hand side allocated"
580 else:
581 print "fix self.__righthandside*=0"
582 self.__righthandside*=0.
583 if self.debug() : print "PDE Debug: Right hand side reset to zero"
584 return self.__righthandside
585
586 #============ some serivice functions =====================================================
587 def getDomain(self):
588 """
589 @brief returns the domain of the PDE
590 """
591 return self.__domain
592
593 def getDim(self):
594 """
595 @brief returns the spatial dimension of the PDE
596 """
597 return self.getDomain().getDim()
598
599 def getNumEquations(self):
600 """
601 @brief returns the number of equations
602 """
603 if self.__numEquations>0:
604 return self.__numEquations
605 else:
606 raise ValueError,"Number of equations is undefined. Please specify argument numEquations."
607
608 def getNumSolutions(self):
609 """
610 @brief returns the number of unknowns
611 """
612 if self.__numSolutions>0:
613 return self.__numSolutions
614 else:
615 raise ValueError,"Number of solution is undefined. Please specify argument numSolutions."
616
617
618 def checkSymmetry(self,verbose=True):
619 """
620 @brief returns if the Operator is symmetric. This is a very expensive operation!!! The symmetry flag is not altered.
621 """
622 verbose=verbose or self.debug()
623 out=True
624 if self.getNumSolutions()!=self.getNumEquations():
625 if verbose: print "non-symmetric PDE because of different number of equations and solutions"
626 out=False
627 else:
628 A=self.getCoefficientOfPDE("A")
629 if not A.isEmpty():
630 tol=util.Lsup(A)*self.TOL
631 if self.getNumSolutions()>1:
632 for i in range(self.getNumEquations()):
633 for j in range(self.getDim()):
634 for k in range(self.getNumSolutions()):
635 for l in range(self.getDim()):
636 if util.Lsup(A[i,j,k,l]-A[k,l,i,j])>tol:
637 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)
638 out=False
639 else:
640 for j in range(self.getDim()):
641 for l in range(self.getDim()):
642 if util.Lsup(A[j,l]-A[l,j])>tol:
643 if verbose: print "non-symmetric PDE because A[%d,%d]!=A[%d,%d]"%(j,l,l,j)
644 out=False
645 B=self.getCoefficientOfPDE("B")
646 C=self.getCoefficientOfPDE("C")
647 if B.isEmpty() and not C.isEmpty():
648 if verbose: print "non-symmetric PDE because B is not present but C is"
649 out=False
650 elif not B.isEmpty() and C.isEmpty():
651 if verbose: print "non-symmetric PDE because C is not present but B is"
652 out=False
653 elif not B.isEmpty() and not C.isEmpty():
654 tol=(util.Lsup(B)+util.Lsup(C))*self.TOL/2.
655 if self.getNumSolutions()>1:
656 for i in range(self.getNumEquations()):
657 for j in range(self.getDim()):
658 for k in range(self.getNumSolutions()):
659 if util.Lsup(B[i,j,k]-C[k,i,j])>tol:
660 if verbose: print "non-symmetric PDE because B[%d,%d,%d]!=C[%d,%d,%d]"%(i,j,k,k,i,j)
661 out=False
662 else:
663 for j in range(self.getDim()):
664 if util.Lsup(B[j]-C[j])>tol:
665 if verbose: print "non-symmetric PDE because B[%d]!=C[%d]"%(j,j)
666 out=False
667 if self.getNumSolutions()>1:
668 D=self.getCoefficientOfPDE("D")
669 if not D.isEmpty():
670 tol=util.Lsup(D)*self.TOL
671 for i in range(self.getNumEquations()):
672 for k in range(self.getNumSolutions()):
673 if util.Lsup(D[i,k]-D[k,i])>tol:
674 if verbose: print "non-symmetric PDE because D[%d,%d]!=D[%d,%d]"%(i,k,k,i)
675 out=False
676
677 return out
678
679 def getFlux(self,u):
680 """
681 @brief returns the flux J_ij for a given u
682
683 J_ij=A_{ijkl}u_{k,l}+B_{ijk}u_k-X_{ij}
684
685 @param u argument of the operator
686
687 """
688 raise SystemError,"getFlux is not implemented yet"
689 return None
690
691 def applyOperator(self,u):
692 """
693 @brief applies the operator of the PDE to a given solution u in weak from
694
695 @param u argument of the operator
696
697 """
698 return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())
699
700 def getResidual(self,u):
701 """
702 @brief return the residual of u in the weak from
703
704 @param u
705 """
706 return self.applyOperator(u)-self.getRightHandSide()
707
708 def _setValue(self,**coefficients):
709 """
710 @brief sets new values to coefficient
711
712 @param coefficients
713 """
714 # check if the coefficients are legal:
715 for i in coefficients.iterkeys():
716 if not self.hasCoefficient(i):
717 raise ValueError,"Attempt to set unknown coefficient %s"%i
718 # if the number of unknowns or equations is still unknown we try to estimate them:
719 if self.__numEquations<1 or self.__numSolutions<1:
720 for i,d in coefficients.iteritems():
721 if hasattr(d,"shape"):
722 s=d.shape
723 elif hasattr(d,"getShape"):
724 s=d.getShape()
725 else:
726 s=numarray.array(d).shape
727 if s!=None:
728 # get number of equations and number of unknowns:
729 res=self.COEFFICIENTS[i].estimateNumEquationsAndNumSolutions(s,self.getDim())
730 if res==None:
731 raise ValueError,"Illegal shape %s of coefficient %s"%(s,i)
732 else:
733 if self.__numEquations<1: self.__numEquations=res[0]
734 if self.__numSolutions<1: self.__numSolutions=res[1]
735 if self.__numEquations<1: raise ValueError,"unidententified number of equations"
736 if self.__numSolutions<1: raise ValueError,"unidententified number of solutions"
737 # now we check the shape of the coefficient if numEquations and numSolutions are set:
738 for i,d in coefficients.iteritems():
739 if d==None:
740 d2=escript.Data()
741 elif isinstance(d,escript.Data):
742 if d.isEmpty():
743 d2=d
744 else:
745 d2=escript.Data(d,self.getFunctionSpaceForCoefficient(i))
746 else:
747 d2=escript.Data(d,self.getFunctionSpaceForCoefficient(i))
748 if not d2.isEmpty():
749 if not self.getShapeOfCoefficient(i)==d2.getShape():
750 raise ValueError,"Expected shape for coefficient %s is %s but actual shape is %s."%(i,self.getShapeOfCoefficient(i),d2.getShape())
751 # overwrite new values:
752 if self.debug(): print "PDE Debug: Coefficient %s has been altered."%i
753 self.COEFFICIENTS[i].setValue(d2)
754 self.alteredCoefficient(i)
755
756 # reset the HomogeneousConstraintFlag:
757 self.__setHomogeneousConstraintFlag()
758 if len(coefficients)>0 and not self.isUsingLumping() and not self.__homogeneous_constraint: self.__rebuildSystem()
759
760 def __setHomogeneousConstraintFlag(self):
761 """
762 @brief checks if the constraints are homogeneous and sets self.__homogeneous_constraint accordingly.
763 """
764 self.__homogeneous_constraint=True
765 q=self.getCoefficientOfPDE("q")
766 r=self.getCoefficientOfPDE("r")
767 if not q.isEmpty() and not r.isEmpty():
768 if (q*r).Lsup()>=1.e-13*r.Lsup(): self.__homogeneous_constraint=False
769 if self.debug():
770 if self.__homogeneous_constraint:
771 print "PDE Debug: Constraints are homogeneous."
772 else:
773 print "PDE Debug: Constraints are inhomogeneous."
774
775
776 # ==== rebuild switches =====================================================================
777 def __rebuildSolution(self,deep=False):
778 """
779 @brief indicates the PDE has to be reolved if the solution is requested
780 """
781 if self.__solution_isValid and self.debug() : print "PDE Debug: PDE has to be resolved."
782 self.__solution_isValid=False
783 if deep: self.__solution=escript.Data()
784
785
786 def __rebuildOperator(self,deep=False):
787 """
788 @brief indicates the operator has to be rebuilt next time it is used
789 """
790 if self.__operator_isValid and self.debug() : print "PDE Debug: Operator has to be rebuilt."
791 self.__rebuildSolution(deep)
792 self.__operator_isValid=False
793 if deep: self.__operator=escript.Operator()
794
795 def __rebuildRightHandSide(self,deep=False):
796 """
797 @brief indicates the right hand side has to be rebuild next time it is used
798 """
799 if self.__righthandside_isValid and self.debug() : print "PDE Debug: Right hand side has to be rebuilt."
800 self.__rebuildSolution(deep)
801 self.__righthandside_isValid=False
802 if deep: self.__righthandside=escript.Data()
803
804 def __rebuildSystem(self,deep=False):
805 """
806 @brief annonced that all coefficient name has been changed
807 """
808 self.__rebuildSolution(deep)
809 self.__rebuildOperator(deep)
810 self.__rebuildRightHandSide(deep)
811
812 def __checkMatrixType(self):
813 """
814 @brief reassess the matrix type and, if needed, initiates an operator rebuild
815 """
816 new_matrix_type=self.getDomain().getSystemMatrixTypeId(self.getSolverMethod(),self.isSymmetric())
817 if not new_matrix_type==self.__matrix_type:
818 if self.debug() : print "PDE Debug: Matrix type is now %d."%new_matrix_type
819 self.__matrix_type=new_matrix_type
820 self.__rebuildOperator(deep=True)
821
822 #============ assembling =======================================================
823 def __copyConstraint(self):
824 """
825 @brief copies the constrint condition into u
826 """
827 if not self.__righthandside.isEmpty():
828 q=self.getCoefficientOfPDE("q")
829 r=self.getCoefficientOfPDE("r")
830 if not q.isEmpty():
831 if r.isEmpty():
832 r2=escript.Data(0,self.__righthandside.getShape(),self.__righthandside.getFunctionSpace())
833 else:
834 r2=escript.Data(r,self.__righthandside.getFunctionSpace())
835 self.__righthandside.copyWithMask(r2,escript.Data(q,self.__righthandside.getFunctionSpace()))
836
837 def __applyConstraint(self):
838 """
839 @brief applies the constraints defined by q and r to the system
840 """
841 q=self.getCoefficientOfPDE("q")
842 r=self.getCoefficientOfPDE("r")
843 if not q.isEmpty() and not self.__operator.isEmpty():
844 # q is the row and column mask to indicate where constraints are set:
845 row_q=escript.Data(q,self.getFunctionSpaceForEquation())
846 col_q=escript.Data(q,self.getFunctionSpaceForSolution())
847 u=self.__makeNewSolution()
848 if r.isEmpty():
849 r_s=self.__makeNewSolution()
850 else:
851 r_s=escript.Data(r,self.getFunctionSpaceForSolution())
852 u.copyWithMask(r_s,col_q)
853 if self.isUsingLumping():
854 self.__operator.copyWithMask(escript.Data(1,q.getShape(),self.getFunctionSpaceForEquation()),row_q)
855 else:
856 if not self.__righthandside.isEmpty(): self.__righthandside-=self.__operator*u
857 self.__operator.nullifyRowsAndCols(row_q,col_q,1.)
858
859 def getSystem(self):
860 """
861 @brief return the operator and right hand side of the PDE
862 """
863 if not self.__operator_isValid or not self.__righthandside_isValid:
864 if self.isUsingLumping():
865 if not self.__operator_isValid:
866 if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():
867 raise TypeError,"Lumped matrix requires same order for equations and unknowns"
868 if not self.getCoefficientOfPDE("A").isEmpty():
869 raise Warning,"Lumped matrix does not allow coefficient A"
870 if not self.getCoefficientOfPDE("B").isEmpty():
871 raise Warning,"Lumped matrix does not allow coefficient B"
872 if not self.getCoefficientOfPDE("C").isEmpty():
873 raise Warning,"Lumped matrix does not allow coefficient C"
874 if self.debug() : print "PDE Debug: New lumped operator is built."
875 mat=self.__makeNewOperator()
876 self.getDomain().addPDEToSystem(mat,escript.Data(), \
877 self.getCoefficientOfPDE("A"), \
878 self.getCoefficientOfPDE("B"), \
879 self.getCoefficientOfPDE("C"), \
880 self.getCoefficientOfPDE("D"), \
881 escript.Data(), \
882 escript.Data(), \
883 self.getCoefficientOfPDE("d"), \
884 escript.Data(),\
885 self.getCoefficientOfPDE("d_contact"), \
886 escript.Data())
887 self.__operator=mat*escript.Data(1,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)
888 self.__applyConstraint()
889 self.__operator_isValid=True
890 if not self.__righthandside_isValid:
891 if self.debug() : print "PDE Debug: New right hand side is built."
892 self.getDomain().addPDEToRHS(self.__getFreshRightHandSide(), \
893 self.getCoefficientOfPDE("X"), \
894 self.getCoefficientOfPDE("Y"),\
895 self.getCoefficientOfPDE("y"),\
896 self.getCoefficientOfPDE("y_contact"))
897 self.__copyConstraint()
898 self.__righthandside_isValid=True
899 else:
900 if not self.__operator_isValid and not self.__righthandside_isValid:
901 if self.debug() : print "PDE Debug: New system is built."
902 self.getDomain().addPDEToSystem(self.__getFreshOperator(),self.__getFreshRightHandSide(), \
903 self.getCoefficientOfPDE("A"), \
904 self.getCoefficientOfPDE("B"), \
905 self.getCoefficientOfPDE("C"), \
906 self.getCoefficientOfPDE("D"), \
907 self.getCoefficientOfPDE("X"), \
908 self.getCoefficientOfPDE("Y"), \
909 self.getCoefficientOfPDE("d"), \
910 self.getCoefficientOfPDE("y"), \
911 self.getCoefficientOfPDE("d_contact"), \
912 self.getCoefficientOfPDE("y_contact"))
913 self.__applyConstraint()
914 self.__copyConstraint()
915 self.__operator_isValid=True
916 self.__righthandside_isValid=True
917 elif not self.__righthandside_isValid:
918 if self.debug() : print "PDE Debug: New right hand side is built."
919 self.getDomain().addPDEToRHS(self.__getFreshRightHandSide(), \
920 self.getCoefficientOfPDE("X"), \
921 self.getCoefficientOfPDE("Y"),\
922 self.getCoefficientOfPDE("y"),\
923 self.getCoefficientOfPDE("y_contact"))
924 self.__copyConstraint()
925 self.__righthandside_isValid=True
926 elif not self.__operator_isValid:
927 if self.debug() : print "PDE Debug: New operator is built."
928 self.getDomain().addPDEToSystem(self.__getFreshOperator(),escript.Data(), \
929 self.getCoefficientOfPDE("A"), \
930 self.getCoefficientOfPDE("B"), \
931 self.getCoefficientOfPDE("C"), \
932 self.getCoefficientOfPDE("D"), \
933 escript.Data(), \
934 escript.Data(), \
935 self.getCoefficientOfPDE("d"), \
936 escript.Data(),\
937 self.getCoefficientOfPDE("d_contact"), \
938 escript.Data())
939 self.__applyConstraint()
940 self.__operator_isValid=True
941 return (self.__operator,self.__righthandside)
942 def getOperator(self):
943 """
944 @brief returns the operator of the PDE
945 """
946 return self.getSystem()[0]
947
948 def getRightHandSide(self):
949 """
950 @brief returns the right hand side of the PDE
951 """
952 return self.getSystem()[1]
953
954 def solve(self,**options):
955 """
956 @brief solve the PDE
957
958 @param options
959 """
960 mat,f=self.getSystem()
961 if self.isUsingLumping():
962 out=f/mat
963 else:
964 options[util.TOLERANCE_KEY]=self.getTolerance()
965 options[util.METHOD_KEY]=self.getSolverMethod()
966 options[util.SYMMETRY_KEY]=self.isSymmetric()
967 if self.debug() : print "PDE Debug: solver options: ",options
968 out=mat.solve(f,options)
969 return out
970
971 def getSolution(self,**options):
972 """
973 @brief returns the solution of the PDE
974
975 @param options
976 """
977 if not self.__solution_isValid:
978 if self.debug() : print "PDE Debug: PDE is resolved."
979 self.__solution=self.solve(**options)
980 self.__solution_isValid=True
981 return self.__solution
982
983
984
985 def ELMAN_RAMAGE(P): return (P-1.).wherePositive()*0.5*(1.-1./(P+1.e-15))
986 def SIMPLIFIED_BROOK_HUGHES(P):
987 c=(P-3.).whereNegative()
988 return P/6.*c+1./2.*(1.-c)
989 def HALF(P): return escript.Scalar(0.5,P.getFunctionSpace())
990
991
992 class AdvectivePDE(LinearPDE):
993 """
994 @brief Class to handel a linear PDE domineated by advective terms:
995
996 class to define a linear PDE of the form
997
998 -(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
999
1000 with boundary conditons:
1001
1002 n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_{ik}u_k = - n_j*X_{ij} + y_i
1003
1004 and contact conditions
1005
1006 n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_contact_{ik}[u_k] = - n_j*X_{ij} + y_contact_i
1007
1008 and constraints:
1009
1010 u_i=r_i where q_i>0
1011
1012 """
1013 def __init__(self,domain,numEquations=0,numSolutions=0,xi=ELMAN_RAMAGE):
1014 LinearPDE.__init__(self,domain,numEquations,numSolutions)
1015 self.__xi=xi
1016 self.__Xi=escript.Data()
1017
1018 def __calculateXi(self,peclet_factor,Z,h):
1019 Z_max=util.Lsup(Z)
1020 if Z_max>0.:
1021 return h*self.__xi(Z*peclet_factor)/(Z+Z_max*self.TOL)
1022 else:
1023 return 0.
1024
1025 def setValue(self,**args):
1026 if "A" in args.keys() or "B" in args.keys() or "C" in args.keys(): self.__Xi=escript.Data()
1027 self._setValue(**args)
1028
1029 def getXi(self):
1030 if self.__Xi.isEmpty():
1031 B=self.getCoefficient("B")
1032 C=self.getCoefficient("C")
1033 A=self.getCoefficient("A")
1034 h=self.getDomain().getSize()
1035 self.__Xi=escript.Scalar(0.,self.getFunctionSpaceForCoefficient("A"))
1036 if not C.isEmpty() or not B.isEmpty():
1037 if not C.isEmpty() and not B.isEmpty():
1038 Z2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))
1039 if self.getNumEquations()>1:
1040 if self.getNumSolutions()>1:
1041 for i in range(self.getNumEquations()):
1042 for k in range(self.getNumSolutions()):
1043 for l in range(self.getDim()): Z2+=(C[i,k,l]-B[i,l,k])**2
1044 else:
1045 for i in range(self.getNumEquations()):
1046 for l in range(self.getDim()): Z2+=(C[i,l]-B[i,l])**2
1047 else:
1048 if self.getNumSolutions()>1:
1049 for k in range(self.getNumSolutions()):
1050 for l in range(self.getDim()): Z2+=(C[k,l]-B[l,k])**2
1051 else:
1052 for l in range(self.getDim()): Z2+=(C[l]-B[l])**2
1053 length_of_Z=util.sqrt(Z2)
1054 elif C.isEmpty():
1055 length_of_Z=util.length(B)
1056 else:
1057 length_of_Z=util.length(C)
1058
1059 Z_max=util.Lsup(length_of_Z)
1060 if Z_max>0.:
1061 length_of_A=util.length(A)
1062 A_max=util.Lsup(length_of_A)
1063 if A_max>0:
1064 inv_A=1./(length_of_A+A_max*self.TOL)
1065 else:
1066 inv_A=1./self.TOL
1067 peclet_number=length_of_Z*h/2*inv_A
1068 xi=self.__xi(peclet_number)
1069 self.__Xi=h*xi/(length_of_Z+Z_max*self.TOL)
1070 print "@ preclet number = %e"%util.Lsup(peclet_number),util.Lsup(xi),util.Lsup(length_of_Z)
1071 return self.__Xi
1072
1073
1074 def getCoefficientOfPDE(self,name):
1075 """
1076 @brief return the value of the coefficient name of the general PDE
1077 @param name
1078 """
1079 if not self.getNumEquations() == self.getNumSolutions():
1080 raise ValueError,"AdvectivePDE expects the number of solution componets and the number of equations to be equal."
1081
1082 if name == "A" :
1083 A=self.getCoefficient("A")
1084 B=self.getCoefficient("B")
1085 C=self.getCoefficient("C")
1086 if B.isEmpty() and C.isEmpty():
1087 Aout=A
1088 else:
1089 if A.isEmpty():
1090 Aout=self.createNewCoefficient("A")
1091 else:
1092 Aout=A[:]
1093 Xi=self.getXi()
1094 if self.getNumEquations()>1:
1095 for i in range(self.getNumEquations()):
1096 for j in range(self.getDim()):
1097 for k in range(self.getNumSolutions()):
1098 for l in range(self.getDim()):
1099 if not C.isEmpty() and not B.isEmpty():
1100 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])
1101 elif C.isEmpty():
1102 for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*B[p,j,i]*B[p,l,k]
1103 else:
1104 for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*C[p,i,j]*C[p,k,l]
1105 else:
1106 for j in range(self.getDim()):
1107 for l in range(self.getDim()):
1108 if not C.isEmpty() and not B.isEmpty():
1109 Aout[j,l]+=Xi*(C[j]-B[j])*(C[l]-B[l])
1110 elif C.isEmpty():
1111 Aout[j,l]+=Xi*B[j]*B[l]
1112 else:
1113 Aout[j,l]+=Xi*C[j]*C[l]
1114 return Aout
1115 elif name == "B" :
1116 B=self.getCoefficient("B")
1117 C=self.getCoefficient("C")
1118 D=self.getCoefficient("D")
1119 if C.isEmpty() or D.isEmpty():
1120 Bout=B
1121 else:
1122 Xi=self.getXi()
1123 if B.isEmpty():
1124 Bout=self.createNewCoefficient("B")
1125 else:
1126 Bout=B[:]
1127 if self.getNumEquations()>1:
1128 for k in range(self.getNumSolutions()):
1129 for p in range(self.getNumEquations()):
1130 tmp=Xi*D[p,k]
1131 for i in range(self.getNumEquations()):
1132 for j in range(self.getDim()):
1133 Bout[i,j,k]+=tmp*C[p,i,j]
1134 else:
1135 tmp=Xi*D
1136 for j in range(self.getDim()): Bout[j]+=tmp*C[j]
1137 return Bout
1138 elif name == "C" :
1139 B=self.getCoefficient("B")
1140 C=self.getCoefficient("C")
1141 D=self.getCoefficient("D")
1142 if B.isEmpty() or D.isEmpty():
1143 Cout=C
1144 else:
1145 Xi=self.getXi()
1146 if C.isEmpty():
1147 Cout=self.createNewCoefficient("C")
1148 else:
1149 Cout=C[:]
1150 if self.getNumEquations()>1:
1151 for k in range(self.getNumSolutions()):
1152 for p in range(self.getNumEquations()):
1153 tmp=Xi*D[p,k]
1154 for i in range(self.getNumEquations()):
1155 for l in range(self.getDim()):
1156 Cout[i,k,l]+=tmp*B[p,l,i]
1157 else:
1158 tmp=Xi*D
1159 for j in range(self.getDim()): Cout[j]+=tmp*B[j]
1160 return Cout
1161 elif name == "D" :
1162 return self.getCoefficient("D")
1163 elif name == "X" :
1164 X=self.getCoefficient("X")
1165 Y=self.getCoefficient("Y")
1166 B=self.getCoefficient("B")
1167 C=self.getCoefficient("C")
1168 if Y.isEmpty() or (B.isEmpty() and C.isEmpty()):
1169 Xout=X
1170 else:
1171 if X.isEmpty():
1172 Xout=self.createNewCoefficient("X")
1173 else:
1174 Xout=X[:]
1175 Xi=self.getXi()
1176 if self.getNumEquations()>1:
1177 for p in range(self.getNumEquations()):
1178 tmp=Xi*Y[p]
1179 for i in range(self.getNumEquations()):
1180 for j in range(self.getDim()):
1181 if not C.isEmpty() and not B.isEmpty():
1182 Xout[i,j]+=tmp*(C[p,i,j]-B[p,j,i])
1183 elif C.isEmpty():
1184 Xout[i,j]-=tmp*B[p,j,i]
1185 else:
1186 Xout[i,j]+=tmp*C[p,i,j]
1187 else:
1188 tmp=Xi*Y
1189 for j in range(self.getDim()):
1190 if not C.isEmpty() and not B.isEmpty():
1191 Xout[j]+=tmp*(C[j]-B[j])
1192 elif C.isEmpty():
1193 Xout[j]-=tmp*B[j]
1194 else:
1195 Xout[j]+=tmp*C[j]
1196 return Xout
1197 elif name == "Y" :
1198 return self.getCoefficient("Y")
1199 elif name == "d" :
1200 return self.getCoefficient("d")
1201 elif name == "y" :
1202 return self.getCoefficient("y")
1203 elif name == "d_contact" :
1204 return self.getCoefficient("d_contact")
1205 elif name == "y_contact" :
1206 return self.getCoefficient("y_contact")
1207 elif name == "r" :
1208 return self.getCoefficient("r")
1209 elif name == "q" :
1210 return self.getCoefficient("q")
1211 else:
1212 raise SystemError,"unknown PDE coefficient %s",name
1213
1214
1215 class Poisson(LinearPDE):
1216 """
1217 @brief Class to define a Poisson equstion problem:
1218
1219 class to define a linear PDE of the form
1220
1221 -u_{,jj} = f
1222
1223 with boundary conditons:
1224
1225 n_j*u_{,j} = 0
1226
1227 and constraints:
1228
1229 u=0 where q>0
1230
1231 """
1232
1233 def __init__(self,domain,f=escript.Data(),q=escript.Data()):
1234 LinearPDE.__init__(self,domain,1,1)
1235 self.COEFFICIENTS={
1236 "f" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1237 "q" : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.BOTH)}
1238 self.setSymmetryOn()
1239 self.setValue(f,q)
1240
1241 def setValue(self,f=escript.Data(),q=escript.Data()):
1242 self._setValue(f=f,q=q)
1243
1244 def getCoefficientOfPDE(self,name):
1245 """
1246 @brief return the value of the coefficient name of the general PDE
1247 @param name
1248 """
1249 if name == "A" :
1250 return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))
1251 elif name == "B" :
1252 return escript.Data()
1253 elif name == "C" :
1254 return escript.Data()
1255 elif name == "D" :
1256 return escript.Data()
1257 elif name == "X" :
1258 return escript.Data()
1259 elif name == "Y" :
1260 return self.getCoefficient("f")
1261 elif name == "d" :
1262 return escript.Data()
1263 elif name == "y" :
1264 return escript.Data()
1265 elif name == "d_contact" :
1266 return escript.Data()
1267 elif name == "y_contact" :
1268 return escript.Data()
1269 elif name == "r" :
1270 return escript.Data()
1271 elif name == "q" :
1272 return self.getCoefficient("q")
1273 else:
1274 raise SystemError,"unknown PDE coefficient %s",name

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26