/[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 304 by gross, Fri Dec 2 06:04:06 2005 UTC revision 425 by gross, Tue Jan 10 04:10:39 2006 UTC
# Line 17  the PDE solver library defined through t Line 17  the PDE solver library defined through t
17  The general interface is provided through the L{LinearPDE} class. The  The general interface is provided through the L{LinearPDE} class. The
18  L{AdvectivePDE} which is derived from the L{LinearPDE} class  L{AdvectivePDE} which is derived from the L{LinearPDE} class
19  provides an interface to PDE dominated by its advective terms. The L{Poisson},  provides an interface to PDE dominated by its advective terms. The L{Poisson},
20  L{Helmholtz}, L{LameEquation}  L{Helmholtz}, L{LameEquation}, L{AdvectionDiffusion}
21  classs which are also derived form the L{LinearPDE} class should be used  classs which are also derived form the L{LinearPDE} class should be used
22  to define of solve these sepecial PDEs.  to define of solve these sepecial PDEs.
23    
# Line 449  class LinearPDE(object): Line 449  class LinearPDE(object):
449     __METHOD_KEY="method"     __METHOD_KEY="method"
450     __SYMMETRY_KEY="symmetric"     __SYMMETRY_KEY="symmetric"
451     __TOLERANCE_KEY="tolerance"     __TOLERANCE_KEY="tolerance"
452       __PRECONDITIONER_KEY="preconditioner"
453    
454    
455     def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):     def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):
# Line 498  class LinearPDE(object): Line 499  class LinearPDE(object):
499       self.__tolerance=1.e-8       self.__tolerance=1.e-8
500       self.__solver_method=self.DEFAULT       self.__solver_method=self.DEFAULT
501       self.__solver_package=self.DEFAULT       self.__solver_package=self.DEFAULT
502         self.__preconditioner=self.DEFAULT
503       self.__matrix_type=self.__domain.getSystemMatrixTypeId(self.DEFAULT,self.DEFAULT,False)       self.__matrix_type=self.__domain.getSystemMatrixTypeId(self.DEFAULT,self.DEFAULT,False)
504       self.__sym=False       self.__sym=False
505    
# Line 772  class LinearPDE(object): Line 774  class LinearPDE(object):
774         @type verbose: C{bool}         @type verbose: C{bool}
775         @keyword reordering: reordering scheme to be used during elimination. Allowed values are         @keyword reordering: reordering scheme to be used during elimination. Allowed values are
776                              L{NO_REORDERING}, L{MINIMUM_FILL_IN}, L{NESTED_DISSECTION}                              L{NO_REORDERING}, L{MINIMUM_FILL_IN}, L{NESTED_DISSECTION}
        @keyword preconditioner: preconditioner method to be used. Allowed values are  
                                 L{SSOR}, L{ILU0}, L{ILUT}, L{JACOBI}  
777         @keyword iter_max: maximum number of iteration steps allowed.         @keyword iter_max: maximum number of iteration steps allowed.
778         @keyword drop_tolerance: threshold for drupping in L{ILUT}         @keyword drop_tolerance: threshold for drupping in L{ILUT}
779         @keyword drop_storage: maximum of allowed memory in L{ILUT}         @keyword drop_storage: maximum of allowed memory in L{ILUT}
# Line 786  class LinearPDE(object): Line 786  class LinearPDE(object):
786               self.__solution=self.copyConstraint(f*mat)               self.__solution=self.copyConstraint(f*mat)
787            else:            else:
788               options[self.__TOLERANCE_KEY]=self.getTolerance()               options[self.__TOLERANCE_KEY]=self.getTolerance()
789               options[self.__METHOD_KEY]=self.getSolverMethod()               options[self.__METHOD_KEY]=self.getSolverMethod()[0]
790                 options[self.__PRECONDITIONER_KEY]=self.getSolverMethod()[1]
791               options[self.__PACKAGE_KEY]=self.getSolverPackage()               options[self.__PACKAGE_KEY]=self.getSolverPackage()
792               options[self.__SYMMETRY_KEY]=self.isSymmetric()               options[self.__SYMMETRY_KEY]=self.isSymmetric()
793               self.trace("PDE is resolved.")               self.trace("PDE is resolved.")
# Line 815  class LinearPDE(object): Line 816  class LinearPDE(object):
816     # =============================================================================     # =============================================================================
817     #   solver settings:     #   solver settings:
818     # =============================================================================     # =============================================================================
819     def setSolverMethod(self,solver=None):     def setSolverMethod(self,solver=None,preconditioner=None):
820         """         """
821         sets a new solver         sets a new solver
822    
823         @param solver: sets a new solver method.         @param solver: sets a new solver method.
824         @type solver: one of L{DEFAULT}, L{ITERATIVE} L{DIRECT}, L{CHOLEVSKY}, L{PCG}, L{CR}, L{CGS}, L{BICGSTAB}, L{SSOR}, L{GMRES}, L{PRES20}, L{LUMPING}.         @type solver: one of L{DEFAULT}, L{ITERATIVE} L{DIRECT}, L{CHOLEVSKY}, L{PCG}, L{CR}, L{CGS}, L{BICGSTAB}, L{SSOR}, L{GMRES}, L{PRES20}, L{LUMPING}.
825           @param preconditioner: sets a new solver method.
826           @type solver: one of L{DEFAULT}, L{JACOBI} L{ILU0}, L{ILUT},L{SSOR}
827         """         """
828         if solver==None: solve=self.DEFAULT         if solver==None: solve=self.DEFAULT
829         if not solver==self.getSolverMethod():         if preconditioner==None: preconditioner=self.DEFAULT
830           if not (solver,preconditioner)==self.getSolverMethod():
831             self.__solver_method=solver             self.__solver_method=solver
832               self.__preconditioner=preconditioner
833             self.__checkMatrixType()             self.__checkMatrixType()
834             self.trace("New solver is %s"%self.getSolverMethodName())             self.trace("New solver is %s"%self.getSolverMethodName())
835    
# Line 838  class LinearPDE(object): Line 843  class LinearPDE(object):
843    
844         m=self.getSolverMethod()         m=self.getSolverMethod()
845         p=self.getSolverPackage()         p=self.getSolverPackage()
846         if m==self.DEFAULT: method="DEFAULT"         method=""
847         elif m==self.DIRECT: method= "DIRECT"         if m[0]==self.DEFAULT: method="DEFAULT"
848         elif m==self.ITERATIVE: method= "ITERATIVE"         elif m[0]==self.DIRECT: method= "DIRECT"
849         elif m==self.CHOLEVSKY: method= "CHOLEVSKY"         elif m[0]==self.ITERATIVE: method= "ITERATIVE"
850         elif m==self.PCG: method= "PCG"         elif m[0]==self.CHOLEVSKY: method= "CHOLEVSKY"
851         elif m==self.CR: method= "CR"         elif m[0]==self.PCG: method= "PCG"
852         elif m==self.CGS: method= "CGS"         elif m[0]==self.CR: method= "CR"
853         elif m==self.BICGSTAB: method= "BICGSTAB"         elif m[0]==self.CGS: method= "CGS"
854         elif m==self.SSOR: method= "SSOR"         elif m[0]==self.BICGSTAB: method= "BICGSTAB"
855         elif m==self.GMRES: method= "GMRES"         elif m[0]==self.SSOR: method= "SSOR"
856         elif m==self.PRES20: method= "PRES20"         elif m[0]==self.GMRES: method= "GMRES"
857         elif m==self.LUMPING: method= "LUMPING"         elif m[0]==self.PRES20: method= "PRES20"
858         else : method="unknown"         elif m[0]==self.LUMPING: method= "LUMPING"
859           if m[1]==self.DEFAULT: method+="+DEFAULT"
860           elif m[1]==self.JACOBI: method+= "+JACOBI"
861           elif m[1]==self.ILU0: method+= "+ILU0"
862           elif m[1]==self.ILUT: method+= "+ILUT"
863           elif m[1]==self.SSOR: method+= "+SSOR"
864         if p==self.DEFAULT: package="DEFAULT"         if p==self.DEFAULT: package="DEFAULT"
865         elif p==self.PASO: package= "PASO"         elif p==self.PASO: package= "PASO"
866         elif p==self.MKL: package= "MKL"         elif p==self.MKL: package= "MKL"
# Line 867  class LinearPDE(object): Line 877  class LinearPDE(object):
877         @return: the solver method currently be used.         @return: the solver method currently be used.
878         @rtype: C{int}         @rtype: C{int}
879         """         """
880         return self.__solver_method         return self.__solver_method,self.__preconditioner
881    
882     def setSolverPackage(self,package=None):     def setSolverPackage(self,package=None):
883         """         """
# Line 898  class LinearPDE(object): Line 908  class LinearPDE(object):
908        @return: True is lumping is currently used a solver method.        @return: True is lumping is currently used a solver method.
909        @rtype: C{bool}        @rtype: C{bool}
910        """        """
911        return self.getSolverMethod()==self.LUMPING        return self.getSolverMethod()[0]==self.LUMPING
912    
913     def setTolerance(self,tol=1.e-8):     def setTolerance(self,tol=1.e-8):
914         """         """
# Line 1091  class LinearPDE(object): Line 1101  class LinearPDE(object):
1101       """       """
1102       reassess the matrix type and, if a new matrix is needed, resets the system.       reassess the matrix type and, if a new matrix is needed, resets the system.
1103       """       """
1104       new_matrix_type=self.getDomain().getSystemMatrixTypeId(self.getSolverMethod(),self.getSolverPackage(),self.isSymmetric())       new_matrix_type=self.getDomain().getSystemMatrixTypeId(self.getSolverMethod()[0],self.getSolverPackage(),self.isSymmetric())
1105       if not new_matrix_type==self.__matrix_type:       if not new_matrix_type==self.__matrix_type:
1106           self.trace("Matrix type is now %d."%new_matrix_type)           self.trace("Matrix type is now %d."%new_matrix_type)
1107           self.__matrix_type=new_matrix_type           self.__matrix_type=new_matrix_type
# Line 1604  class Poisson(LinearPDE): Line 1614  class Poisson(LinearPDE):
1614    
1615     """     """
1616    
1617     def __init__(self,domain,f=escript.Data(),q=escript.Data(),debug=False):     def __init__(self,domain,debug=False):
1618       """       """
1619       initializes a new Poisson equation       initializes a new Poisson equation
1620    
# Line 1644  class Poisson(LinearPDE): Line 1654  class Poisson(LinearPDE):
1654       @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.
1655       """       """
1656       if name == "A" :       if name == "A" :
1657           return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))           return escript.Data(util.kronecker(self.getDim()),escript.Function(self.getDomain()))
1658       elif name == "B" :       elif name == "B" :
1659           return escript.Data()           return escript.Data()
1660       elif name == "C" :       elif name == "C" :
# Line 2018  class AdvectivePDE(LinearPDE): Line 2028  class AdvectivePDE(LinearPDE):
2028       """       """
2029       return escript.Scalar(0.5,P.getFunctionSpace())       return escript.Scalar(0.5,P.getFunctionSpace())
2030    
2031     def __calculateXi(self,peclet_factor,Z,h):     def __calculateXi(self,peclet_factor,flux,h):
2032         Z_max=util.Lsup(Z)         flux=util.Lsup(flux)
2033         if Z_max>0.:         if flux_max>0.:
2034            return h*self.__xi(Z*peclet_factor)/(Z+Z_max*self.__TOL)            return h*self.__xi(flux*peclet_factor)/(flux+flux_max*self.__TOL)
2035         else:         else:
2036            return 0.            return 0.
2037    
# Line 2034  class AdvectivePDE(LinearPDE): Line 2044  class AdvectivePDE(LinearPDE):
2044           self.__Xi=escript.Scalar(0.,self.getFunctionSpaceForCoefficient("A"))           self.__Xi=escript.Scalar(0.,self.getFunctionSpaceForCoefficient("A"))
2045           if not C.isEmpty() or not B.isEmpty():           if not C.isEmpty() or not B.isEmpty():
2046              if not C.isEmpty() and not B.isEmpty():              if not C.isEmpty() and not B.isEmpty():
2047                  Z2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))                  flux2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))
2048                  if self.getNumEquations()>1:                  if self.getNumEquations()>1:
2049                     if self.getNumSolutions()>1:                     if self.getNumSolutions()>1:
2050                        for i in range(self.getNumEquations()):                        for i in range(self.getNumEquations()):
2051                           for k in range(self.getNumSolutions()):                           for k in range(self.getNumSolutions()):
2052                              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
2053                          # flux=C-util.reorderComponents(B,[0,2,1])
2054                     else:                     else:
2055                        for i in range(self.getNumEquations()):                        for i in range(self.getNumEquations()):
2056                           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
2057                          # flux=C-B
2058                  else:                  else:
2059                     if self.getNumSolutions()>1:                     if self.getNumSolutions()>1:
2060                        for k in range(self.getNumSolutions()):                        for k in range(self.getNumSolutions()):
2061                           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
2062                          # flux=C-util.reorderComponents(B,[1,0])
2063                     else:                     else:
2064                        for l in range(self.getDim()): Z2+=(C[l]-B[l])**2                        for l in range(self.getDim()): flux2+=(C[l]-B[l])**2
2065                  length_of_Z=util.sqrt(Z2)                        #flux=C-B
2066                    length_of_flux=util.sqrt(flux2)
2067              elif C.isEmpty():              elif C.isEmpty():
2068                length_of_Z=util.length(B)                length_of_flux=util.length(B)
2069                  #flux=B
2070              else:              else:
2071                length_of_Z=util.length(C)                length_of_flux=util.length(C)
2072                  #flux=C
2073    
2074              Z_max=util.Lsup(length_of_Z)              #length_of_flux=util.length(flux)
2075              if Z_max>0.:              flux_max=util.Lsup(length_of_flux)
2076                if flux_max>0.:
2077                   # length_of_A=util.inner(flux,util.tensormutiply(A,flux))
2078                 length_of_A=util.length(A)                 length_of_A=util.length(A)
2079                 A_max=util.Lsup(length_of_A)                 A_max=util.Lsup(length_of_A)
2080                 if A_max>0:                 if A_max>0:
2081                      inv_A=1./(length_of_A+A_max*self.__TOL)                      inv_A=1./(length_of_A+A_max*self.__TOL)
2082                 else:                 else:
2083                      inv_A=1./self.__TOL                      inv_A=1./self.__TOL
2084                 peclet_number=length_of_Z*h/2*inv_A                 peclet_number=length_of_flux*h/2*inv_A
2085                 xi=self.__xi(peclet_number)                 xi=self.__xi(peclet_number)
2086                 self.__Xi=h*xi/(length_of_Z+Z_max*self.__TOL)                 self.__Xi=h*xi/(length_of_flux+flux_max*self.__TOL)
2087                 self.trace("preclet number = %e"%util.Lsup(peclet_number))                 self.trace("preclet number = %e"%util.Lsup(peclet_number))
2088        return self.__Xi        return self.__Xi
2089    
# Line 2103  class AdvectivePDE(LinearPDE): Line 2121  class AdvectivePDE(LinearPDE):
2121                        for k in range(self.getNumSolutions()):                        for k in range(self.getNumSolutions()):
2122                           for l in range(self.getDim()):                           for l in range(self.getDim()):
2123                              if not C.isEmpty() and not B.isEmpty():                              if not C.isEmpty() and not B.isEmpty():
2124                                   # tmp=C-util.reorderComponents(B,[0,2,1])
2125                                   # Aout=Aout+Xi*util.generalTensorProduct(util.reorder(tmp,[1,2,0]),tmp,offset=1)
2126                                 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])
2127                              elif C.isEmpty():                              elif C.isEmpty():
2128                                 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]
2129                                   # Aout=Aout+Xi*util.generalTensorProduct(util.reorder(B,[2,1,0]),util.reorder(B,[0,2,1]),offset=1)
2130                              else:                              else:
2131                                 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]
2132                                   # Aout=Aout+Xi*util.generalTensorProduct(util.reorder(C,[1,2,0]),C,offset=1)
2133              else:              else:
2134                  for j in range(self.getDim()):                  for j in range(self.getDim()):
2135                     for l in range(self.getDim()):                     for l in range(self.getDim()):
# Line 2117  class AdvectivePDE(LinearPDE): Line 2139  class AdvectivePDE(LinearPDE):
2139                            Aout[j,l]+=Xi*B[j]*B[l]                            Aout[j,l]+=Xi*B[j]*B[l]
2140                        else:                        else:
2141                            Aout[j,l]+=Xi*C[j]*C[l]                            Aout[j,l]+=Xi*C[j]*C[l]
2142                     # if not C.isEmpty() and not B.isEmpty():
2143                     #    tmp=C-B
2144                     #    Aout=Aout+Xi*util.outer(tmp,tmp)
2145                     # elif C.isEmpty():
2146                     #    Aout=Aout+Xi*util.outer(B,B)
2147                     # else:
2148                     # Aout=Aout+Xi*util.outer(C,C)
2149           return Aout           return Aout
2150       elif name == "B" :       elif name == "B" :
2151           B=self.getCoefficient("B")           B=self.getCoefficient("B")
# Line 2137  class AdvectivePDE(LinearPDE): Line 2166  class AdvectivePDE(LinearPDE):
2166                       for i in range(self.getNumEquations()):                       for i in range(self.getNumEquations()):
2167                          for j in range(self.getDim()):                          for j in range(self.getDim()):
2168                             Bout[i,j,k]+=tmp*C[p,i,j]                             Bout[i,j,k]+=tmp*C[p,i,j]
2169                               # Bout=Bout+Xi*util.generalTensorProduct(util.reorder(C,[1,2,0]),D,offset=1)
2170              else:              else:
2171                 tmp=Xi*D                 tmp=Xi*D
2172                 for j in range(self.getDim()): Bout[j]+=tmp*C[j]                 for j in range(self.getDim()): Bout[j]+=tmp*C[j]
2173                   # Bout=Bout+Xi*D*C
2174           return Bout           return Bout
2175       elif name == "C" :       elif name == "C" :
2176           B=self.getCoefficient("B")           B=self.getCoefficient("B")
# Line 2160  class AdvectivePDE(LinearPDE): Line 2191  class AdvectivePDE(LinearPDE):
2191                        for i in range(self.getNumEquations()):                        for i in range(self.getNumEquations()):
2192                          for l in range(self.getDim()):                          for l in range(self.getDim()):
2193                                   Cout[i,k,l]+=tmp*B[p,l,i]                                   Cout[i,k,l]+=tmp*B[p,l,i]
2194                                     # Cout=Cout+Xi*B[p,l,i]*D[p,k]
2195              else:              else:
2196                 tmp=Xi*D                 tmp=Xi*D
2197                 for j in range(self.getDim()): Cout[j]+=tmp*B[j]                 for j in range(self.getDim()): Cout[j]+=tmp*B[j]
2198                   # Cout=Cout+tmp*D*B
2199           return Cout           return Cout
2200       elif name == "D" :       elif name == "D" :
2201           return self.getCoefficient("D")           return self.getCoefficient("D")
# Line 2186  class AdvectivePDE(LinearPDE): Line 2219  class AdvectivePDE(LinearPDE):
2219                         for j in range(self.getDim()):                         for j in range(self.getDim()):
2220                            if not C.isEmpty() and not B.isEmpty():                            if not C.isEmpty() and not B.isEmpty():
2221                               Xout[i,j]+=tmp*(C[p,i,j]-B[p,j,i])                               Xout[i,j]+=tmp*(C[p,i,j]-B[p,j,i])
2222                                 # Xout=X_out+Xi*util.inner(Y,C-util.reorderComponents(B,[0,2,1]),offset=1)
2223                            elif C.isEmpty():                            elif C.isEmpty():
2224                               Xout[i,j]-=tmp*B[p,j,i]                               Xout[i,j]-=tmp*B[p,j,i]
2225                                 # Xout=X_out-Xi*util.inner(Y,util.reorderComponents(B,[0,2,1]),offset=1)
2226                            else:                            else:
2227                               Xout[i,j]+=tmp*C[p,i,j]                               Xout[i,j]+=tmp*C[p,i,j]
2228                                 # Xout=X_out+Xi*util.inner(Y,C,offset=1)
2229              else:              else:
2230                   tmp=Xi*Y                   tmp=Xi*Y
2231                   for j in range(self.getDim()):                   for j in range(self.getDim()):
2232                      if not C.isEmpty() and not B.isEmpty():                      if not C.isEmpty() and not B.isEmpty():
2233                         Xout[j]+=tmp*(C[j]-B[j])                         Xout[j]+=tmp*(C[j]-B[j])
2234                           # Xout=Xout+Xi*Y*(C-B)
2235                      elif C.isEmpty():                      elif C.isEmpty():
2236                         Xout[j]-=tmp*B[j]                         Xout[j]-=tmp*B[j]
2237                           # Xout=Xout-Xi*Y*B
2238                      else:                      else:
2239                         Xout[j]+=tmp*C[j]                         Xout[j]+=tmp*C[j]
2240                           # Xout=Xout+Xi*Y*C
2241           return Xout           return Xout
2242       elif name == "Y" :       elif name == "Y" :
2243           return self.getCoefficient("Y")           return self.getCoefficient("Y")
# Line 2213  class AdvectivePDE(LinearPDE): Line 2252  class AdvectivePDE(LinearPDE):
2252       elif name == "r" :       elif name == "r" :
2253           return self.getCoefficient("r")           return self.getCoefficient("r")
2254       elif name == "q" :       elif name == "q" :
2255             return self.getCoefficient("q")
2256         else:
2257            raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
2258    
2259    class AdvectionDiffusion(LinearPDE):
2260       """
2261       Class to define PDE equation of the unisotropic advection-diffusion problem, which is genear L{LinearPDE} of the form
2262    
2263       M{S{omega}*u + inner(v,grad(u))- grad(matrixmult(k_bar,grad(u))[j])[j] = f}
2264    
2265       with natural boundary conditons
2266    
2267       M{n[j]*matrixmult(k,grad(u))[j] = g- S{alpha}u }
2268    
2269       and constraints:
2270    
2271       M{u=r} where M{q>0}
2272    
2273       and
2274    
2275       M{k_bar[i,j]=k[i,j]+upwind[i]*upwind[j]}
2276    
2277       """
2278    
2279       def __init__(self,domain,debug=False):
2280         """
2281         initializes a new Poisson equation
2282    
2283         @param domain: domain of the PDE
2284         @type domain: L{Domain<escript.Domain>}
2285         @param debug: if True debug informations are printed.
2286    
2287         """
2288         super(AdvectionDiffusion, self).__init__(domain,1,1,debug)
2289         self.COEFFICIENTS={"omega": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
2290                            "k": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_DIM,PDECoefficient.BY_DIM),PDECoefficient.OPERATOR),
2291                            "f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
2292                            "v": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_DIM,),PDECoefficient.OPERATOR),
2293                            "upwind": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_DIM,),PDECoefficient.OPERATOR),
2294                            "alpha": PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
2295                            "g": PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
2296                            "r": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH),
2297                            "q": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}
2298    
2299       def setValue(self,**coefficients):
2300         """
2301         sets new values to coefficients
2302    
2303         @param coefficients: new values assigned to coefficients
2304         @keyword omega: value for coefficient M{S{omega}}
2305         @type omega: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
2306         @keyword k: value for coefficient M{k}
2307         @type k: any type that can be casted to L{Tensor<escript.Tensor>} object on L{Function<escript.Function>}.
2308         @keyword v: value for coefficient M{v}
2309         @type v: any type that can be casted to L{Vector<escript.Vector>} object on L{Function<escript.Function>}.
2310         @keyword upwind: value for upwind term M{upwind}
2311         @type upwind: any type that can be casted to L{Vector<escript.Vector>} object on L{Function<escript.Function>}.
2312         @keyword f: value for right hand side M{f}
2313         @type f: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
2314         @keyword alpha: value for right hand side M{S{alpha}}
2315         @type alpha: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
2316         @keyword g: value for right hand side M{g}
2317         @type g: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
2318         @keyword r: prescribed values M{r} for the solution in constraints.
2319         @type r: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2320                   depending of reduced order is used for the representation of the equation.
2321         @keyword q: mask for location of constraints
2322         @type q: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2323                   depending of reduced order is used for the representation of the equation.
2324         @raise IllegalCoefficient: if an unknown coefficient keyword is used.
2325         """
2326         super(AdvectionDiffusion, self).setValue(**coefficients)
2327    
2328       def getCoefficientOfGeneralPDE(self,name):
2329         """
2330         return the value of the coefficient name of the general PDE
2331    
2332         @param name: name of the coefficient requested.
2333         @type name: C{string}
2334         @return: the value of the coefficient  name
2335         @rtype: L{Data<escript.Data>}
2336         @raise IllegalCoefficient: if name is not one of coefficients
2337                      "A", M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact}, M{r} or M{q}.
2338         @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
2339         """
2340         if name == "A" :
2341             return self.getCoefficient("k")+util.outer(self.getCoefficient("upwind"),self.getCoefficient("upwind"))
2342         elif name == "B" :
2343             return escript.Data()
2344         elif name == "C" :
2345             return self.getCoefficient("v")
2346         elif name == "D" :
2347             return self.getCoefficient("omega")
2348         elif name == "X" :
2349             return escript.Data()
2350         elif name == "Y" :
2351             return self.getCoefficient("f")
2352         elif name == "d" :
2353             return self.getCoefficient("alpha")
2354         elif name == "y" :
2355             return self.getCoefficient("g")
2356         elif name == "d_contact" :
2357             return escript.Data()
2358         elif name == "y_contact" :
2359             return escript.Data()
2360         elif name == "r" :
2361             return self.getCoefficient("r")
2362         elif name == "q" :
2363           return self.getCoefficient("q")           return self.getCoefficient("q")
2364       else:       else:
2365          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name

Legend:
Removed from v.304  
changed lines
  Added in v.425

  ViewVC Help
Powered by ViewVC 1.1.26