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

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

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

revision 327 by gross, Fri Dec 2 06:04:06 2005 UTC revision 328 by gross, Wed Dec 7 04:41:53 2005 UTC
# Line 1604  class Poisson(LinearPDE): Line 1604  class Poisson(LinearPDE):
1604    
1605     """     """
1606    
1607     def __init__(self,domain,f=escript.Data(),q=escript.Data(),debug=False):     def __init__(self,domain,debug=False):
1608       """       """
1609       initializes a new Poisson equation       initializes a new Poisson equation
1610    
# Line 1644  class Poisson(LinearPDE): Line 1644  class Poisson(LinearPDE):
1644       @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.       @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
1645       """       """
1646       if name == "A" :       if name == "A" :
1647           return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))           return escript.Data(util.kronecker(self.getDim()),escript.Function(self.getDomain()))
1648       elif name == "B" :       elif name == "B" :
1649           return escript.Data()           return escript.Data()
1650       elif name == "C" :       elif name == "C" :
# Line 2018  class AdvectivePDE(LinearPDE): Line 2018  class AdvectivePDE(LinearPDE):
2018       """       """
2019       return escript.Scalar(0.5,P.getFunctionSpace())       return escript.Scalar(0.5,P.getFunctionSpace())
2020    
2021     def __calculateXi(self,peclet_factor,Z,h):     def __calculateXi(self,peclet_factor,flux,h):
2022         Z_max=util.Lsup(Z)         flux=util.Lsup(flux)
2023         if Z_max>0.:         if flux_max>0.:
2024            return h*self.__xi(Z*peclet_factor)/(Z+Z_max*self.__TOL)            return h*self.__xi(flux*peclet_factor)/(flux+flux_max*self.__TOL)
2025         else:         else:
2026            return 0.            return 0.
2027    
# Line 2034  class AdvectivePDE(LinearPDE): Line 2034  class AdvectivePDE(LinearPDE):
2034           self.__Xi=escript.Scalar(0.,self.getFunctionSpaceForCoefficient("A"))           self.__Xi=escript.Scalar(0.,self.getFunctionSpaceForCoefficient("A"))
2035           if not C.isEmpty() or not B.isEmpty():           if not C.isEmpty() or not B.isEmpty():
2036              if not C.isEmpty() and not B.isEmpty():              if not C.isEmpty() and not B.isEmpty():
2037                  Z2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))                  flux2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))
2038                  if self.getNumEquations()>1:                  if self.getNumEquations()>1:
2039                     if self.getNumSolutions()>1:                     if self.getNumSolutions()>1:
2040                        for i in range(self.getNumEquations()):                        for i in range(self.getNumEquations()):
2041                           for k in range(self.getNumSolutions()):                           for k in range(self.getNumSolutions()):
2042                              for l in range(self.getDim()): Z2+=(C[i,k,l]-B[i,l,k])**2                              for l in range(self.getDim()): flux2+=(C[i,k,l]-B[i,l,k])**2
2043                          # flux=C-util.reorderComponents(B,[0,2,1])
2044                     else:                     else:
2045                        for i in range(self.getNumEquations()):                        for i in range(self.getNumEquations()):
2046                           for l in range(self.getDim()): Z2+=(C[i,l]-B[i,l])**2                           for l in range(self.getDim()): flux2+=(C[i,l]-B[i,l])**2
2047                          # flux=C-B
2048                  else:                  else:
2049                     if self.getNumSolutions()>1:                     if self.getNumSolutions()>1:
2050                        for k in range(self.getNumSolutions()):                        for k in range(self.getNumSolutions()):
2051                           for l in range(self.getDim()): Z2+=(C[k,l]-B[l,k])**2                           for l in range(self.getDim()): flux2+=(C[k,l]-B[l,k])**2
2052                          # flux=C-util.reorderComponents(B,[1,0])
2053                     else:                     else:
2054                        for l in range(self.getDim()): Z2+=(C[l]-B[l])**2                        for l in range(self.getDim()): flux2+=(C[l]-B[l])**2
2055                  length_of_Z=util.sqrt(Z2)                        #flux=C-B
2056                    length_of_flux=util.sqrt(flux2)
2057              elif C.isEmpty():              elif C.isEmpty():
2058                length_of_Z=util.length(B)                length_of_flux=util.length(B)
2059                  #flux=B
2060              else:              else:
2061                length_of_Z=util.length(C)                length_of_flux=util.length(C)
2062                  #flux=C
2063    
2064              Z_max=util.Lsup(length_of_Z)              #length_of_flux=util.length(flux)
2065              if Z_max>0.:              flux_max=util.Lsup(length_of_flux)
2066                if flux_max>0.:
2067                   # length_of_A=util.inner(flux,util.tensormutiply(A,flux))
2068                 length_of_A=util.length(A)                 length_of_A=util.length(A)
2069                 A_max=util.Lsup(length_of_A)                 A_max=util.Lsup(length_of_A)
2070                 if A_max>0:                 if A_max>0:
2071                      inv_A=1./(length_of_A+A_max*self.__TOL)                      inv_A=1./(length_of_A+A_max*self.__TOL)
2072                 else:                 else:
2073                      inv_A=1./self.__TOL                      inv_A=1./self.__TOL
2074                 peclet_number=length_of_Z*h/2*inv_A                 peclet_number=length_of_flux*h/2*inv_A
2075                 xi=self.__xi(peclet_number)                 xi=self.__xi(peclet_number)
2076                 self.__Xi=h*xi/(length_of_Z+Z_max*self.__TOL)                 self.__Xi=h*xi/(length_of_flux+flux_max*self.__TOL)
2077                 self.trace("preclet number = %e"%util.Lsup(peclet_number))                 self.trace("preclet number = %e"%util.Lsup(peclet_number))
2078        return self.__Xi        return self.__Xi
2079    
# Line 2103  class AdvectivePDE(LinearPDE): Line 2111  class AdvectivePDE(LinearPDE):
2111                        for k in range(self.getNumSolutions()):                        for k in range(self.getNumSolutions()):
2112                           for l in range(self.getDim()):                           for l in range(self.getDim()):
2113                              if not C.isEmpty() and not B.isEmpty():                              if not C.isEmpty() and not B.isEmpty():
2114                                   # tmp=C-util.reorderComponents(B,[0,2,1])
2115                                   # Aout=Aout+Xi*util.generalTensorProduct(util.reorder(tmp,[1,2,0]),tmp,offset=1)
2116                                 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])                                 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])
2117                              elif C.isEmpty():                              elif C.isEmpty():
2118                                 for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*B[p,j,i]*B[p,l,k]                                 for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*B[p,j,i]*B[p,l,k]
2119                                   # Aout=Aout+Xi*util.generalTensorProduct(util.reorder(B,[2,1,0]),util.reorder(B,[0,2,1]),offset=1)
2120                              else:                              else:
2121                                 for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*C[p,i,j]*C[p,k,l]                                 for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*C[p,i,j]*C[p,k,l]
2122                                   # Aout=Aout+Xi*util.generalTensorProduct(util.reorder(C,[1,2,0]),C,offset=1)
2123              else:              else:
2124                  for j in range(self.getDim()):                  for j in range(self.getDim()):
2125                     for l in range(self.getDim()):                     for l in range(self.getDim()):
# Line 2117  class AdvectivePDE(LinearPDE): Line 2129  class AdvectivePDE(LinearPDE):
2129                            Aout[j,l]+=Xi*B[j]*B[l]                            Aout[j,l]+=Xi*B[j]*B[l]
2130                        else:                        else:
2131                            Aout[j,l]+=Xi*C[j]*C[l]                            Aout[j,l]+=Xi*C[j]*C[l]
2132                     # if not C.isEmpty() and not B.isEmpty():
2133                     #    tmp=C-B
2134                     #    Aout=Aout+Xi*util.outer(tmp,tmp)
2135                     # elif C.isEmpty():
2136                     #    Aout=Aout+Xi*util.outer(B,B)
2137                     # else:
2138                     # Aout=Aout+Xi*util.outer(C,C)
2139           return Aout           return Aout
2140       elif name == "B" :       elif name == "B" :
2141           B=self.getCoefficient("B")           B=self.getCoefficient("B")
# Line 2137  class AdvectivePDE(LinearPDE): Line 2156  class AdvectivePDE(LinearPDE):
2156                       for i in range(self.getNumEquations()):                       for i in range(self.getNumEquations()):
2157                          for j in range(self.getDim()):                          for j in range(self.getDim()):
2158                             Bout[i,j,k]+=tmp*C[p,i,j]                             Bout[i,j,k]+=tmp*C[p,i,j]
2159                               # Bout=Bout+Xi*util.generalTensorProduct(util.reorder(C,[1,2,0]),D,offset=1)
2160              else:              else:
2161                 tmp=Xi*D                 tmp=Xi*D
2162                 for j in range(self.getDim()): Bout[j]+=tmp*C[j]                 for j in range(self.getDim()): Bout[j]+=tmp*C[j]
2163                   # Bout=Bout+Xi*D*C
2164           return Bout           return Bout
2165       elif name == "C" :       elif name == "C" :
2166           B=self.getCoefficient("B")           B=self.getCoefficient("B")
# Line 2160  class AdvectivePDE(LinearPDE): Line 2181  class AdvectivePDE(LinearPDE):
2181                        for i in range(self.getNumEquations()):                        for i in range(self.getNumEquations()):
2182                          for l in range(self.getDim()):                          for l in range(self.getDim()):
2183                                   Cout[i,k,l]+=tmp*B[p,l,i]                                   Cout[i,k,l]+=tmp*B[p,l,i]
2184                                     # Cout=Cout+Xi*B[p,l,i]*D[p,k]
2185              else:              else:
2186                 tmp=Xi*D                 tmp=Xi*D
2187                 for j in range(self.getDim()): Cout[j]+=tmp*B[j]                 for j in range(self.getDim()): Cout[j]+=tmp*B[j]
2188                   # Cout=Cout+tmp*D*B
2189           return Cout           return Cout
2190       elif name == "D" :       elif name == "D" :
2191           return self.getCoefficient("D")           return self.getCoefficient("D")
# Line 2186  class AdvectivePDE(LinearPDE): Line 2209  class AdvectivePDE(LinearPDE):
2209                         for j in range(self.getDim()):                         for j in range(self.getDim()):
2210                            if not C.isEmpty() and not B.isEmpty():                            if not C.isEmpty() and not B.isEmpty():
2211                               Xout[i,j]+=tmp*(C[p,i,j]-B[p,j,i])                               Xout[i,j]+=tmp*(C[p,i,j]-B[p,j,i])
2212                                 # Xout=X_out+Xi*util.inner(Y,C-util.reorderComponents(B,[0,2,1]),offset=1)
2213                            elif C.isEmpty():                            elif C.isEmpty():
2214                               Xout[i,j]-=tmp*B[p,j,i]                               Xout[i,j]-=tmp*B[p,j,i]
2215                                 # Xout=X_out-Xi*util.inner(Y,util.reorderComponents(B,[0,2,1]),offset=1)
2216                            else:                            else:
2217                               Xout[i,j]+=tmp*C[p,i,j]                               Xout[i,j]+=tmp*C[p,i,j]
2218                                 # Xout=X_out+Xi*util.inner(Y,C,offset=1)
2219              else:              else:
2220                   tmp=Xi*Y                   tmp=Xi*Y
2221                   for j in range(self.getDim()):                   for j in range(self.getDim()):
2222                      if not C.isEmpty() and not B.isEmpty():                      if not C.isEmpty() and not B.isEmpty():
2223                         Xout[j]+=tmp*(C[j]-B[j])                         Xout[j]+=tmp*(C[j]-B[j])
2224                           # Xout=Xout+Xi*Y*(C-B)
2225                      elif C.isEmpty():                      elif C.isEmpty():
2226                         Xout[j]-=tmp*B[j]                         Xout[j]-=tmp*B[j]
2227                           # Xout=Xout-Xi*Y*B
2228                      else:                      else:
2229                         Xout[j]+=tmp*C[j]                         Xout[j]+=tmp*C[j]
2230                           # Xout=Xout+Xi*Y*C
2231           return Xout           return Xout
2232       elif name == "Y" :       elif name == "Y" :
2233           return self.getCoefficient("Y")           return self.getCoefficient("Y")
Legend:
Removed from v.327  
changed lines
  Added in v.328
  ViewVC Help
Powered by ViewVC 1.1.26