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

revision 108 by jgs, Thu Jan 27 06:21:59 2005 UTC revision 114 by jgs, Fri Mar 4 07:12:37 2005 UTC
# Line 172  class LinearPDE: Line 172  class LinearPDE:
172     SSOR=util.SSOR     SSOR=util.SSOR
173     GMRES=util.GMRES     GMRES=util.GMRES
174     PRES20=util.PRES20     PRES20=util.PRES20
175       ILU0=util.ILU0
176       JACOBI=util.JACOBI
177
178     def __init__(self,domain,numEquations=0,numSolutions=0):     def __init__(self,domain,numEquations=0,numSolutions=0):
179       """       """
# Line 224  class LinearPDE: Line 226  class LinearPDE:
226       @param name       @param name
227       """       """
228       return escript.Data(shape = getShapeOfCoefficient(name), \       return escript.Data(shape = getShapeOfCoefficient(name), \
229                           what = getFunctionSpaceOfCoefficient(name))                           what = getFunctionSpaceForCoefficient(name))
230
231     def __del__(self):     def __del__(self):
232       pass       pass
# Line 656  class LinearPDE: Line 658  class LinearPDE:
658                 for i in range(self.getNumEquations()):                 for i in range(self.getNumEquations()):
659                     for j in range(self.getDim()):                     for j in range(self.getDim()):
660                        for k in range(self.getNumSolutions()):                        for k in range(self.getNumSolutions()):
661                           if util.Lsup(B[i,j,k]-C[i,j,k])>tol:                           if util.Lsup(B[i,j,k]-C[k,i,j])>tol:
662                                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)
663                                out=False                                out=False
664              else:              else:
665                 for j in range(self.getDim()):                 for j in range(self.getDim()):
# Line 980  class LinearPDE: Line 982  class LinearPDE:
982             self.__solution_isValid=True             self.__solution_isValid=True
983         return self.__solution         return self.__solution
984
985
986
987    def ELMAN_RAMAGE(P): return (P-1.).wherePositive()*0.5*(1.-1./(P+1.e-15))
988    def SIMPLIFIED_BROOK_HUGHES(P):
989             c=(P-3.).whereNegative()
990             return P/6.*c+1./2.*(1.-c)
991    def HALF(P): return escript.Scalar(0.5,P.getFunctionSpace())
992
993
995     """     """
996     @brief Class to handel a linear PDE domineated by advective terms:     @brief Class to handel a linear PDE domineated by advective terms:
1011
1012           u_i=r_i where q_i>0           u_i=r_i where q_i>0
1013
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
1014     """     """
1015     def __getXi(self,alpha):     def __init__(self,domain,numEquations=0,numSolutions=0,xi=ELMAN_RAMAGE):
1016           c=alpha-3.        LinearPDE.__init__(self,domain,numEquations,numSolutions)
1017           return c*c.whereNegative()/3.+1.        self.__xi=xi
1018          self.__Xi=escript.Data()
1019     def __getUpdateVector(self,V,hover2,alphaByU):
1020       v=util.length(V)     def __calculateXi(self,peclet_factor,Z,h):
1021       v_max=util.Lsup(v)         Z_max=util.Lsup(Z)
1022       if v_max>0:         if Z_max>0.:
1023           V/=v+v_max*self.TOL            return h*self.__xi(Z*peclet_factor)/(Z+Z_max*self.TOL)
1024           alpha=alphaByU*v         else:
1025           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
1026
1027     def __getAlphaByU(self,A,hover2):     def setValue(self,**args):
1028        a=util.length(A)         if "A" in args.keys()   or "B" in args.keys() or "C" in args.keys(): self.__Xi=escript.Data()
1029        a_max=util.Lsup(a)         self._setValue(**args)
1030        if a_max>0:
1031           return hover2/(a+a_max*self.TOL)     def getXi(self):
1032        else:        if self.__Xi.isEmpty():
1033           return 1./self.TOL           B=self.getCoefficient("B")
1034             C=self.getCoefficient("C")
1035             A=self.getCoefficient("A")
1036             h=self.getDomain().getSize()
1037             self.__Xi=escript.Scalar(0.,self.getFunctionSpaceForCoefficient("A"))
1038             if not C.isEmpty() or not B.isEmpty():
1039                if not C.isEmpty() and not B.isEmpty():
1040                    Z2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))
1041                    if self.getNumEquations()>1:
1042                       if self.getNumSolutions()>1:
1043                          for i in range(self.getNumEquations()):
1044                             for k in range(self.getNumSolutions()):
1045                                for l in range(self.getDim()): Z2+=(C[i,k,l]-B[i,l,k])**2
1046                       else:
1047                          for i in range(self.getNumEquations()):
1048                             for l in range(self.getDim()): Z2+=(C[i,l]-B[i,l])**2
1049                    else:
1050                       if self.getNumSolutions()>1:
1051                          for k in range(self.getNumSolutions()):
1052                             for l in range(self.getDim()): Z2+=(C[k,l]-B[l,k])**2
1053                       else:
1054                          for l in range(self.getDim()): Z2+=(C[l]-B[l])**2
1055                    length_of_Z=util.sqrt(Z2)
1056                elif C.isEmpty():
1057                  length_of_Z=util.length(B)
1058                else:
1059                  length_of_Z=util.length(C)
1060
1061                Z_max=util.Lsup(length_of_Z)
1062                if Z_max>0.:
1063                   length_of_A=util.length(A)
1064                   A_max=util.Lsup(length_of_A)
1065                   if A_max>0:
1066                        inv_A=1./(length_of_A+A_max*self.TOL)
1067                   else:
1068                        inv_A=1./self.TOL
1069                   peclet_number=length_of_Z*h/2*inv_A
1070                   xi=self.__xi(peclet_number)
1071                   self.__Xi=h*xi/(length_of_Z+Z_max*self.TOL)
1072                   print "@ preclet number = %e"%util.Lsup(peclet_number),util.Lsup(xi),util.Lsup(length_of_Z)
1073          return self.__Xi
1074
1075
1076     def getCoefficientOfPDE(self,name):     def getCoefficientOfPDE(self,name):
1077       """       """
1078       @brief return the value of the coefficient name of the general PDE       @brief return the value of the coefficient name of the general PDE
1079       @param name       @param name
1080       """       """
1081         if not self.getNumEquations() == self.getNumSolutions():
1082              raise ValueError,"AdvectivePDE expects the number of solution componets and the number of equations to be equal."
1083
1084       if name == "A" :       if name == "A" :
1085           A=self.getCoefficient("A")           A=self.getCoefficient("A")
1086           B=self.getCoefficient("B")           B=self.getCoefficient("B")
1087           C=self.getCoefficient("C")           C=self.getCoefficient("C")
1088           if not B.isEmpty() or not C.isEmpty():           if B.isEmpty() and C.isEmpty():
1089               if A.isEmpty():              Aout=A
1090                   A=self.createNewCoefficient("A")           else:
1091               else:              if A.isEmpty():
1092                   A=A[:]                 Aout=self.createNewCoefficient("A")
1093               hover2=self.getDomain().getSize()/2.              else:
1094               if self.getNumEquations()>1:                 Aout=A[:]
1095                  if self.getNumSolutions()>1:              Xi=self.getXi()
1096                     for i in range(self.getNumEquations()):              if self.getNumEquations()>1:
1097                    for i in range(self.getNumEquations()):
1098                       for j in range(self.getDim()):
1099                        for k in range(self.getNumSolutions()):                        for k in range(self.getNumSolutions()):
1100                           alphaByU=self.__getAlphaByU(A[i,:,k,:],hover2)                           for l in range(self.getDim()):
1101                           if not B.isEmpty():                              if not C.isEmpty() and not B.isEmpty():
1102                               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])
1103                               for j in range(self.getDim()):                              elif C.isEmpty():
1104                                  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]
1105                                     A[i,j,k,l]+=f*b_sub[j]*b_sub[l]                              else:
1106                           if not C.isEmpty():                                 for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*C[p,i,j]*C[p,k,l]
1107                               c_sub,f=self.__getUpdateVector(C[i,k,:],hover2,alphaByU)              else:
1108                               for j in range(self.getDim()):                  for j in range(self.getDim()):
1109                                  for l in range(self.getDim()):                     for l in range(self.getDim()):
1110                                     A[i,j,k,l]+=f*c_sub[j]*c_sub[l]                        if not C.isEmpty() and not B.isEmpty():
1111                  else:                              Aout[j,l]+=Xi*(C[j]-B[j])*(C[l]-B[l])
1112                     for i in range(self.getNumEquations()):                        elif C.isEmpty():
1113                        alphaByU=self.__getAlphaByU(A[i,:,:],hover2)                            Aout[j,l]+=Xi*B[j]*B[l]
1114                        if not B.isEmpty():                        else:
1115                            b_sub,f=self.__getUpdateVector(B[i,:],hover2,alphaByU)                            Aout[j,l]+=Xi*C[j]*C[l]
1116                            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
1117       elif name == "B" :       elif name == "B" :
1118           return self.getCoefficient("B")           B=self.getCoefficient("B")
1119             C=self.getCoefficient("C")
1120             D=self.getCoefficient("D")
1121             if C.isEmpty() or D.isEmpty():
1122                Bout=B
1123             else:
1124                Xi=self.getXi()
1125                if B.isEmpty():
1126                    Bout=self.createNewCoefficient("B")
1127                else:
1128                    Bout=B[:]
1129                if self.getNumEquations()>1:
1130                   for k in range(self.getNumSolutions()):
1131                      for p in range(self.getNumEquations()):
1132                         tmp=Xi*D[p,k]
1133                         for i in range(self.getNumEquations()):
1134                            for j in range(self.getDim()):
1135                               Bout[i,j,k]+=tmp*C[p,i,j]
1136                else:
1137                   tmp=Xi*D
1138                   for j in range(self.getDim()): Bout[j]+=tmp*C[j]
1139             return Bout
1140       elif name == "C" :       elif name == "C" :
1141           return self.getCoefficient("C")           B=self.getCoefficient("B")
1142             C=self.getCoefficient("C")
1143             D=self.getCoefficient("D")
1144             if B.isEmpty() or D.isEmpty():
1145                Cout=C
1146             else:
1147                Xi=self.getXi()
1148                if C.isEmpty():
1149                    Cout=self.createNewCoefficient("C")
1150                else:
1151                    Cout=C[:]
1152                if self.getNumEquations()>1:
1153                   for k in range(self.getNumSolutions()):
1154                       for p in range(self.getNumEquations()):
1155                          tmp=Xi*D[p,k]
1156                          for i in range(self.getNumEquations()):
1157                            for l in range(self.getDim()):
1158                                     Cout[i,k,l]+=tmp*B[p,l,i]
1159                else:
1160                   tmp=Xi*D
1161                   for j in range(self.getDim()): Cout[j]+=tmp*B[j]
1162             return Cout
1163       elif name == "D" :       elif name == "D" :
1164           return self.getCoefficient("D")           return self.getCoefficient("D")
1165       elif name == "X" :       elif name == "X" :
1166           return self.getCoefficient("X")           X=self.getCoefficient("X")
1167             Y=self.getCoefficient("Y")
1168             B=self.getCoefficient("B")
1169             C=self.getCoefficient("C")
1170             if Y.isEmpty() or (B.isEmpty() and C.isEmpty()):
1171                Xout=X
1172             else:
1173                if X.isEmpty():
1174                    Xout=self.createNewCoefficient("X")
1175                else:
1176                    Xout=X[:]
1177                Xi=self.getXi()
1178                if self.getNumEquations()>1:
1179                     for p in range(self.getNumEquations()):
1180                        tmp=Xi*Y[p]
1181                        for i in range(self.getNumEquations()):
1182                           for j in range(self.getDim()):
1183                              if not C.isEmpty() and not B.isEmpty():
1184                                 Xout[i,j]+=tmp*(C[p,i,j]-B[p,j,i])
1185                              elif C.isEmpty():
1186                                 Xout[i,j]-=tmp*B[p,j,i]
1187                              else:
1188                                 Xout[i,j]+=tmp*C[p,i,j]
1189                else:
1190                     tmp=Xi*Y
1191                     for j in range(self.getDim()):
1192                        if not C.isEmpty() and not B.isEmpty():
1193                           Xout[j]+=tmp*(C[j]-B[j])
1194                        elif C.isEmpty():
1195                           Xout[j]-=tmp*B[j]
1196                        else:
1197                           Xout[j]+=tmp*C[j]
1198             return Xout
1199       elif name == "Y" :       elif name == "Y" :
1200           return self.getCoefficient("Y")           return self.getCoefficient("Y")
1201       elif name == "d" :       elif name == "d" :

Legend:
 Removed from v.108 changed lines Added in v.114