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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 108 by jgs, Thu Jan 27 06:21:59 2005 UTC revision 110 by jgs, Mon Feb 14 04:14:42 2005 UTC
# Line 224  class LinearPDE: Line 224  class LinearPDE:
224       @param name       @param name
225       """       """
226       return escript.Data(shape = getShapeOfCoefficient(name), \       return escript.Data(shape = getShapeOfCoefficient(name), \
227                           what = getFunctionSpaceOfCoefficient(name))                           what = getFunctionSpaceForCoefficient(name))
228    
229     def __del__(self):     def __del__(self):
230       pass       pass
# Line 656  class LinearPDE: Line 656  class LinearPDE:
656                 for i in range(self.getNumEquations()):                 for i in range(self.getNumEquations()):
657                     for j in range(self.getDim()):                     for j in range(self.getDim()):
658                        for k in range(self.getNumSolutions()):                        for k in range(self.getNumSolutions()):
659                           if util.Lsup(B[i,j,k]-C[i,j,k])>tol:                           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,i,j,k)                                if verbose: print "non-symmetric PDE because B[%d,%d,%d]!=C[%d,%d,%d]"%(i,j,k,k,i,j)
661                                out=False                                out=False
662              else:              else:
663                 for j in range(self.getDim()):                 for j in range(self.getDim()):
# Line 980  class LinearPDE: Line 980  class LinearPDE:
980             self.__solution_isValid=True             self.__solution_isValid=True
981         return self.__solution         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):  class AdvectivePDE(LinearPDE):
993     """     """
994     @brief Class to handel a linear PDE domineated by advective terms:     @brief Class to handel a linear PDE domineated by advective terms:
# Line 1000  class AdvectivePDE(LinearPDE): Line 1009  class AdvectivePDE(LinearPDE):
1009    
1010           u_i=r_i where q_i>0           u_i=r_i where q_i>0
1011    
     The PDE is solved by stabilizing the advective terms using SUPG approach:  
   
        A_{ijkl}<-A_{ijkl}+0.5*h*(xi(b_{ik})*B_{ijk}*B_{ilk}/length(B_{i:k})^2)+0.5*h*xi_{c_{ik}}*(C_{ikj}*C_{ikl}/length(C_{ik:})^2)  
   
     where    
   
            b_{ik}=length(B_{i:k})*h/2/length(A_{i:k:})  
            c_{ik}=length(C_{i:k})*h/2/length(A_{i:k:})  
   
                       alpha/3        alpha<3  
            xi(alpha)=          for                  approximating cotanh(alpha)-1/alpha  
                        1             alpha>=3  
1012     """     """
1013     def __getXi(self,alpha):     def __init__(self,domain,numEquations=0,numSolutions=0,xi=ELMAN_RAMAGE):
1014           c=alpha-3.        LinearPDE.__init__(self,domain,numEquations,numSolutions)
1015           return c*c.whereNegative()/3.+1.        self.__xi=xi
1016          self.__Xi=escript.Data()
1017     def __getUpdateVector(self,V,hover2,alphaByU):  
1018       v=util.length(V)     def __calculateXi(self,peclet_factor,Z,h):
1019       v_max=util.Lsup(v)         Z_max=util.Lsup(Z)
1020       if v_max>0:         if Z_max>0.:
1021           V/=v+v_max*self.TOL            return h*self.__xi(Z*peclet_factor)/(Z+Z_max*self.TOL)
1022           alpha=alphaByU*v         else:
1023           A_bar=v*hover2*self.__getXi(alpha)            return 0.
          print "-------------"  
          print "@ max alpha ",util.Lsup(alpha)  
          print "-------------"  
      else:  
          A_bar=1.  
      return V,A_bar  
1024    
1025     def __getAlphaByU(self,A,hover2):     def setValue(self,**args):
1026        a=util.length(A)         if "A" in args.keys()   or "B" in args.keys() or "C" in args.keys(): self.__Xi=escript.Data()
1027        a_max=util.Lsup(a)         self._setValue(**args)
1028        if a_max>0:            
1029           return hover2/(a+a_max*self.TOL)     def getXi(self):
1030        else:        if self.__Xi.isEmpty():
1031           return 1./self.TOL           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):     def getCoefficientOfPDE(self,name):
1075       """       """
1076       @brief return the value of the coefficient name of the general PDE       @brief return the value of the coefficient name of the general PDE
1077       @param name       @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" :       if name == "A" :
1083           A=self.getCoefficient("A")           A=self.getCoefficient("A")
1084           B=self.getCoefficient("B")           B=self.getCoefficient("B")
1085           C=self.getCoefficient("C")           C=self.getCoefficient("C")
1086           if not B.isEmpty() or not C.isEmpty():           if B.isEmpty() and C.isEmpty():
1087               if A.isEmpty():              Aout=A
1088                   A=self.createNewCoefficient("A")           else:
1089               else:              if A.isEmpty():
1090                   A=A[:]                 Aout=self.createNewCoefficient("A")
1091               hover2=self.getDomain().getSize()/2.              else:
1092               if self.getNumEquations()>1:                 Aout=A[:]
1093                  if self.getNumSolutions()>1:              Xi=self.getXi()
1094                     for i in range(self.getNumEquations()):              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()):                        for k in range(self.getNumSolutions()):
1098                           alphaByU=self.__getAlphaByU(A[i,:,k,:],hover2)                           for l in range(self.getDim()):
1099                           if not B.isEmpty():                              if not C.isEmpty() and not B.isEmpty():
1100                               b_sub,f=self.__getUpdateVector(B[i,:,k],hover2,alphaByU)                                 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                               for j in range(self.getDim()):                              elif C.isEmpty():
1102                                  for l in range(self.getDim()):                                 for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*B[p,j,i]*B[p,l,k]
1103                                     A[i,j,k,l]+=f*b_sub[j]*b_sub[l]                              else:
1104                           if not C.isEmpty():                                 for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*C[p,i,j]*C[p,k,l]
1105                               c_sub,f=self.__getUpdateVector(C[i,k,:],hover2,alphaByU)              else:
1106                               for j in range(self.getDim()):                  for j in range(self.getDim()):
1107                                  for l in range(self.getDim()):                     for l in range(self.getDim()):
1108                                     A[i,j,k,l]+=f*c_sub[j]*c_sub[l]                        if not C.isEmpty() and not B.isEmpty():
1109                  else:                              Aout[j,l]+=Xi*(C[j]-B[j])*(C[l]-B[l])
1110                     for i in range(self.getNumEquations()):                        elif C.isEmpty():
1111                        alphaByU=self.__getAlphaByU(A[i,:,:],hover2)                            Aout[j,l]+=Xi*B[j]*B[l]
1112                        if not B.isEmpty():                        else:
1113                            b_sub,f=self.__getUpdateVector(B[i,:],hover2,alphaByU)                            Aout[j,l]+=Xi*C[j]*C[l]
1114                            for j in range(self.getDim()):           return Aout
                              for l in range(self.getDim()):  
                                  A[i,j,l]+=f*b_sub[j]*b_sub[l]  
                       if not C.isEmpty():  
                            c_sub,f=self.__getUpdateVector(C[i,:],hover2,alphaByU)  
                            for j in range(self.getDim()):  
                               for l in range(self.getDim()):  
                                  A[i,j,l]+=f*c_sub[j]*c_sub[l]  
              else:  
                 if self.getNumSolutions()>1:  
                    for k in range(self.getNumSolutions()):  
                       alphaByU=self.__getAlphaByU(A[:,k,:],hover2)  
                       if not B.isEmpty():  
                          b_sub,f=self.__getUpdateVector(B[:,k],hover2,alphaByU)  
                          for j in range(self.getDim()):  
                             for l in range(self.getDim()):  
                                    A[j,k,l]+=f*b_sub[j]*b_sub[l]  
                       if not C.isEmpty():  
                          c_sub,f=self.__getUpdateVector(C[k,:],hover2,alphaByU)  
                          for j in range(self.getDim()):  
                             for l in range(self.getDim()):  
                                A[j,k,l]+=f*c_sub[j]*c_sub[l]  
                 else:    
                    alphaByU=self.__getAlphaByU(A[:,:],hover2)  
                    if not B.isEmpty():  
                        b_sub,f=self.__getUpdateVector(B[:],hover2,alphaByU)  
                        for j in range(self.getDim()):  
                           for l in range(self.getDim()):  
                              A[j,l]+=f*b_sub[j]*b_sub[l]  
                    if not C.isEmpty():  
                       c_sub,f=self.__getUpdateVector(C[:],hover2,alphaByU)  
                       for j in range(self.getDim()):  
                           for l in range(self.getDim()):  
                              A[j,l]+=f*c_sub[j]*c_sub[l]  
          return A  
1115       elif name == "B" :       elif name == "B" :
1116           return self.getCoefficient("B")           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" :       elif name == "C" :
1139           return self.getCoefficient("C")           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" :       elif name == "D" :
1162           return self.getCoefficient("D")           return self.getCoefficient("D")
1163       elif name == "X" :       elif name == "X" :
1164           return self.getCoefficient("X")           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" :       elif name == "Y" :
1198           return self.getCoefficient("Y")           return self.getCoefficient("Y")
1199       elif name == "d" :       elif name == "d" :
Legend:
Removed from v.108  
changed lines
  Added in v.110
  ViewVC Help
Powered by ViewVC 1.1.26