# Annotation of /trunk/esys2/escript/py_src/linearPDEs.py

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

```
 1 jgs 102 # \$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 jgs 108 class PDECoefficient: 27 jgs 102 """ 28 @brief 29 """ 30 jgs 108 # identifier for location of Data objects defining COEFFICIENTS 31 jgs 102 INTERIOR=0 32 BOUNDARY=1 33 CONTACT=2 34 CONTINUOUS=3 35 jgs 108 # identifier in the pattern of COEFFICIENTS: 36 jgs 102 # 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 jgs 108 self.resetValue() 53 jgs 102 54 jgs 108 def resetValue(self): 55 """ 56 @brief resets coefficient value to default 57 """ 58 self.value=escript.Data() 59 60 jgs 102 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 jgs 108 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 jgs 102 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 jgs 108 @brief Class to handel a linear PDE 146 jgs 102 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 jgs 108 TOL=1.e-13 165 jgs 102 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 jgs 104 def __init__(self,domain,numEquations=0,numSolutions=0): 177 jgs 102 """ 178 @brief initializes a new linear PDE. 179 180 @param args 181 """ 182 jgs 108 # 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 jgs 102 197 # initialize attributes 198 self.__debug=None 199 jgs 104 self.__domain=domain 200 self.__numEquations=numEquations 201 self.__numSolutions=numSolutions 202 jgs 102 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 jgs 108 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 jgs 110 what = getFunctionSpaceForCoefficient(name)) 228 jgs 108 229 def __del__(self): 230 pass 231 232 jgs 102 def getCoefficient(self,name): 233 """ 234 jgs 108 @brief return the value of the parameter name 235 jgs 102 236 @param name 237 """ 238 jgs 108 return self.COEFFICIENTS[name].getValue() 239 jgs 102 240 jgs 108 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 jgs 102 """ 250 jgs 108 @brief return true if name is the name of a coefficient 251 jgs 102 252 jgs 108 @param name 253 jgs 102 """ 254 jgs 108 return self.COEFFICIENTS.has_key(name) 255 jgs 102 256 jgs 108 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 jgs 102 """ 270 @brief sets new values to coefficients 271 272 @param coefficients 273 """ 274 jgs 108 self._setValue(**coefficients) 275 jgs 102 276 277 def cleanCoefficients(self): 278 """ 279 @brief resets all coefficients to default values. 280 """ 281 jgs 108 for i in self.COEFFICIENTS.iterkeys(): 282 self.COEFFICIENTS[i].resetValue() 283 jgs 102 284 jgs 108 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 jgs 102 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 jgs 108 return self.COEFFICIENTS[name].buildShape(self.getNumEquations(),self.getNumSolutions(),self.getDomain().getDim()) 299 jgs 102 else: 300 raise ValueError,"Solution coefficient %s requested"%name 301 302 jgs 108 def getFunctionSpaceForCoefficient(self,name): 303 jgs 102 """ 304 @brief return the atoms of the coefficient name 305 306 @param name 307 """ 308 if self.hasCoefficient(name): 309 jgs 108 return self.COEFFICIENTS[name].getFunctionSpace(self.getDomain()) 310 jgs 102 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 jgs 108 if self.COEFFICIENTS[name].isAlteringOperator(): self.__rebuildOperator() 321 if self.COEFFICIENTS[name].isAlteringRightHandSide(): self.__rebuildRightHandSide() 322 jgs 102 else: 323 jgs 108 raise ValueError,"unknown coefficient %s requested"%name 324 jgs 102 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 tol0: 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 jgs 110 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 jgs 108 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 jgs 102 # ==== 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 jgs 108 if deep: self.__solution=escript.Data() 784 jgs 102 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 jgs 108 def __copyConstraint(self): 824 jgs 102 """ 825 @brief copies the constrint condition into u 826 """ 827 jgs 108 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 jgs 102 837 jgs 108 def __applyConstraint(self): 838 jgs 102 """ 839 jgs 108 @brief applies the constraints defined by q and r to the system 840 jgs 102 """ 841 jgs 108 q=self.getCoefficientOfPDE("q") 842 r=self.getCoefficientOfPDE("r") 843 jgs 102 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 jgs 108 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 jgs 102 859 jgs 108 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 jgs 102 def getOperator(self): 943 """ 944 @brief returns the operator of the PDE 945 """ 946 jgs 108 return self.getSystem()[0] 947 jgs 102 948 jgs 108 def getRightHandSide(self): 949 jgs 102 """ 950 @brief returns the right hand side of the PDE 951 """ 952 jgs 108 return self.getSystem()[1] 953 jgs 102 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 jgs 110 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 jgs 108 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 jgs 104 998 jgs 108 -(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 jgs 102 1000 jgs 108 with boundary conditons: 1001 jgs 102 1002 jgs 108 n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_{ik}u_k = - n_j*X_{ij} + y_i 1003 jgs 102 1004 jgs 108 and contact conditions 1005 jgs 102 1006 jgs 108 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 jgs 102 1008 jgs 108 and constraints: 1009 jgs 102 1010 jgs 108 u_i=r_i where q_i>0 1011 jgs 102 1012 jgs 110 """ 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 jgs 102 1018 jgs 110 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 jgs 102 1025 jgs 110 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 jgs 102 1059 jgs 110 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 jgs 102 1074 jgs 108 def getCoefficientOfPDE(self,name): 1075 """ 1076 @brief return the value of the coefficient name of the general PDE 1077 @param name 1078 """ 1079 jgs 110 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 jgs 108 if name == "A" : 1083 A=self.getCoefficient("A") 1084 B=self.getCoefficient("B") 1085 C=self.getCoefficient("C") 1086 jgs 110 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 jgs 108 for k in range(self.getNumSolutions()): 1098 jgs 110 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 jgs 108 elif name == "B" : 1116 jgs 110 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 jgs 108 elif name == "C" : 1139 jgs 110 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 jgs 108 elif name == "D" : 1162 return self.getCoefficient("D") 1163 elif name == "X" : 1164 jgs 110 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 jgs 108 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 jgs 102 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 jgs 108 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 jgs 102 self.setSymmetryOn() 1239 self.setValue(f,q) 1240 1241 def setValue(self,f=escript.Data(),q=escript.Data()): 1242 jgs 108 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