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

Diff of /temp_trunk_copy/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 430 by gross, Wed Jan 11 06:40:50 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 419  class LinearPDE(object): Line 419  class LinearPDE(object):
419     @cvar MKL: Intel's MKL solver library     @cvar MKL: Intel's MKL solver library
420     @cvar UMFPACK: the UMFPACK library     @cvar UMFPACK: the UMFPACK library
421     @cvar ITERATIVE: The default iterative solver     @cvar ITERATIVE: The default iterative solver
422       @cvar AMG: algebraic multi grid
423       @cvar RILU: recursive ILU
424    
425     """     """
426     DEFAULT= 0     DEFAULT= 0
# Line 443  class LinearPDE(object): Line 445  class LinearPDE(object):
445     UMFPACK= 16     UMFPACK= 16
446     ITERATIVE= 20     ITERATIVE= 20
447     PASO= 21     PASO= 21
448       AMG= 22
449       RILU = 23
450    
451     __TOL=1.e-13     __TOL=1.e-13
452     __PACKAGE_KEY="package"     __PACKAGE_KEY="package"
453     __METHOD_KEY="method"     __METHOD_KEY="method"
454     __SYMMETRY_KEY="symmetric"     __SYMMETRY_KEY="symmetric"
455     __TOLERANCE_KEY="tolerance"     __TOLERANCE_KEY="tolerance"
456       __PRECONDITIONER_KEY="preconditioner"
457    
458    
459     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 503  class LinearPDE(object):
503       self.__tolerance=1.e-8       self.__tolerance=1.e-8
504       self.__solver_method=self.DEFAULT       self.__solver_method=self.DEFAULT
505       self.__solver_package=self.DEFAULT       self.__solver_package=self.DEFAULT
506         self.__preconditioner=self.DEFAULT
507       self.__matrix_type=self.__domain.getSystemMatrixTypeId(self.DEFAULT,self.DEFAULT,False)       self.__matrix_type=self.__domain.getSystemMatrixTypeId(self.DEFAULT,self.DEFAULT,False)
508       self.__sym=False       self.__sym=False
509    
# Line 772  class LinearPDE(object): Line 778  class LinearPDE(object):
778         @type verbose: C{bool}         @type verbose: C{bool}
779         @keyword reordering: reordering scheme to be used during elimination. Allowed values are         @keyword reordering: reordering scheme to be used during elimination. Allowed values are
780                              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}  
781         @keyword iter_max: maximum number of iteration steps allowed.         @keyword iter_max: maximum number of iteration steps allowed.
782         @keyword drop_tolerance: threshold for drupping in L{ILUT}         @keyword drop_tolerance: threshold for drupping in L{ILUT}
783         @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 790  class LinearPDE(object):
790               self.__solution=self.copyConstraint(f*mat)               self.__solution=self.copyConstraint(f*mat)
791            else:            else:
792               options[self.__TOLERANCE_KEY]=self.getTolerance()               options[self.__TOLERANCE_KEY]=self.getTolerance()
793               options[self.__METHOD_KEY]=self.getSolverMethod()               options[self.__METHOD_KEY]=self.getSolverMethod()[0]
794                 options[self.__PRECONDITIONER_KEY]=self.getSolverMethod()[1]
795               options[self.__PACKAGE_KEY]=self.getSolverPackage()               options[self.__PACKAGE_KEY]=self.getSolverPackage()
796               options[self.__SYMMETRY_KEY]=self.isSymmetric()               options[self.__SYMMETRY_KEY]=self.isSymmetric()
797               self.trace("PDE is resolved.")               self.trace("PDE is resolved.")
# Line 815  class LinearPDE(object): Line 820  class LinearPDE(object):
820     # =============================================================================     # =============================================================================
821     #   solver settings:     #   solver settings:
822     # =============================================================================     # =============================================================================
823     def setSolverMethod(self,solver=None):     def setSolverMethod(self,solver=None,preconditioner=None):
824         """         """
825         sets a new solver         sets a new solver
826    
827         @param solver: sets a new solver method.         @param solver: sets a new solver method.
828         @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}.
829           @param preconditioner: sets a new solver method.
830           @type solver: one of L{DEFAULT}, L{JACOBI} L{ILU0}, L{ILUT},L{SSOR}
831         """         """
832         if solver==None: solve=self.DEFAULT         if solver==None: solve=self.DEFAULT
833         if not solver==self.getSolverMethod():         if preconditioner==None: preconditioner=self.DEFAULT
834           if not (solver,preconditioner)==self.getSolverMethod():
835             self.__solver_method=solver             self.__solver_method=solver
836               self.__preconditioner=preconditioner
837             self.__checkMatrixType()             self.__checkMatrixType()
838             self.trace("New solver is %s"%self.getSolverMethodName())             self.trace("New solver is %s"%self.getSolverMethodName())
839    
# Line 838  class LinearPDE(object): Line 847  class LinearPDE(object):
847    
848         m=self.getSolverMethod()         m=self.getSolverMethod()
849         p=self.getSolverPackage()         p=self.getSolverPackage()
850         if m==self.DEFAULT: method="DEFAULT"         method=""
851         elif m==self.DIRECT: method= "DIRECT"         if m[0]==self.DEFAULT: method="DEFAULT"
852         elif m==self.ITERATIVE: method= "ITERATIVE"         elif m[0]==self.DIRECT: method= "DIRECT"
853         elif m==self.CHOLEVSKY: method= "CHOLEVSKY"         elif m[0]==self.ITERATIVE: method= "ITERATIVE"
854         elif m==self.PCG: method= "PCG"         elif m[0]==self.CHOLEVSKY: method= "CHOLEVSKY"
855         elif m==self.CR: method= "CR"         elif m[0]==self.PCG: method= "PCG"
856         elif m==self.CGS: method= "CGS"         elif m[0]==self.CR: method= "CR"
857         elif m==self.BICGSTAB: method= "BICGSTAB"         elif m[0]==self.CGS: method= "CGS"
858         elif m==self.SSOR: method= "SSOR"         elif m[0]==self.BICGSTAB: method= "BICGSTAB"
859         elif m==self.GMRES: method= "GMRES"         elif m[0]==self.SSOR: method= "SSOR"
860         elif m==self.PRES20: method= "PRES20"         elif m[0]==self.GMRES: method= "GMRES"
861         elif m==self.LUMPING: method= "LUMPING"         elif m[0]==self.PRES20: method= "PRES20"
862         else : method="unknown"         elif m[0]==self.LUMPING: method= "LUMPING"
863           if m[1]==self.DEFAULT: method+="+DEFAULT"
864           elif m[1]==self.JACOBI: method+= "+JACOBI"
865           elif m[1]==self.ILU0: method+= "+ILU0"
866           elif m[1]==self.ILUT: method+= "+ILUT"
867           elif m[1]==self.SSOR: method+= "+SSOR"
868         if p==self.DEFAULT: package="DEFAULT"         if p==self.DEFAULT: package="DEFAULT"
869         elif p==self.PASO: package= "PASO"         elif p==self.PASO: package= "PASO"
870         elif p==self.MKL: package= "MKL"         elif p==self.MKL: package= "MKL"
# Line 867  class LinearPDE(object): Line 881  class LinearPDE(object):
881         @return: the solver method currently be used.         @return: the solver method currently be used.
882         @rtype: C{int}         @rtype: C{int}
883         """         """
884         return self.__solver_method         return self.__solver_method,self.__preconditioner
885    
886     def setSolverPackage(self,package=None):     def setSolverPackage(self,package=None):
887         """         """
# Line 898  class LinearPDE(object): Line 912  class LinearPDE(object):
912        @return: True is lumping is currently used a solver method.        @return: True is lumping is currently used a solver method.
913        @rtype: C{bool}        @rtype: C{bool}
914        """        """
915        return self.getSolverMethod()==self.LUMPING        return self.getSolverMethod()[0]==self.LUMPING
916    
917     def setTolerance(self,tol=1.e-8):     def setTolerance(self,tol=1.e-8):
918         """         """
# Line 1091  class LinearPDE(object): Line 1105  class LinearPDE(object):
1105       """       """
1106       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.
1107       """       """
1108       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())
1109       if not new_matrix_type==self.__matrix_type:       if not new_matrix_type==self.__matrix_type:
1110           self.trace("Matrix type is now %d."%new_matrix_type)           self.trace("Matrix type is now %d."%new_matrix_type)
1111           self.__matrix_type=new_matrix_type           self.__matrix_type=new_matrix_type
# Line 1604  class Poisson(LinearPDE): Line 1618  class Poisson(LinearPDE):
1618    
1619     """     """
1620    
1621     def __init__(self,domain,f=escript.Data(),q=escript.Data(),debug=False):     def __init__(self,domain,debug=False):
1622       """       """
1623       initializes a new Poisson equation       initializes a new Poisson equation
1624    
# Line 1644  class Poisson(LinearPDE): Line 1658  class Poisson(LinearPDE):
1658       @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.
1659       """       """
1660       if name == "A" :       if name == "A" :
1661           return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))           return escript.Data(util.kronecker(self.getDim()),escript.Function(self.getDomain()))
1662       elif name == "B" :       elif name == "B" :
1663           return escript.Data()           return escript.Data()
1664       elif name == "C" :       elif name == "C" :
# Line 2018  class AdvectivePDE(LinearPDE): Line 2032  class AdvectivePDE(LinearPDE):
2032       """       """
2033       return escript.Scalar(0.5,P.getFunctionSpace())       return escript.Scalar(0.5,P.getFunctionSpace())
2034    
2035     def __calculateXi(self,peclet_factor,Z,h):     def __calculateXi(self,peclet_factor,flux,h):
2036         Z_max=util.Lsup(Z)         flux=util.Lsup(flux)
2037         if Z_max>0.:         if flux_max>0.:
2038            return h*self.__xi(Z*peclet_factor)/(Z+Z_max*self.__TOL)            return h*self.__xi(flux*peclet_factor)/(flux+flux_max*self.__TOL)
2039         else:         else:
2040            return 0.            return 0.
2041    
# Line 2034  class AdvectivePDE(LinearPDE): Line 2048  class AdvectivePDE(LinearPDE):
2048           self.__Xi=escript.Scalar(0.,self.getFunctionSpaceForCoefficient("A"))           self.__Xi=escript.Scalar(0.,self.getFunctionSpaceForCoefficient("A"))
2049           if not C.isEmpty() or not B.isEmpty():           if not C.isEmpty() or not B.isEmpty():
2050              if not C.isEmpty() and not B.isEmpty():              if not C.isEmpty() and not B.isEmpty():
2051                  Z2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))                  flux2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))
2052                  if self.getNumEquations()>1:                  if self.getNumEquations()>1:
2053                     if self.getNumSolutions()>1:                     if self.getNumSolutions()>1:
2054                        for i in range(self.getNumEquations()):                        for i in range(self.getNumEquations()):
2055                           for k in range(self.getNumSolutions()):                           for k in range(self.getNumSolutions()):
2056                              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
2057                          # flux=C-util.reorderComponents(B,[0,2,1])
2058                     else:                     else:
2059                        for i in range(self.getNumEquations()):                        for i in range(self.getNumEquations()):
2060                           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
2061                          # flux=C-B
2062                  else:                  else:
2063                     if self.getNumSolutions()>1:                     if self.getNumSolutions()>1:
2064                        for k in range(self.getNumSolutions()):                        for k in range(self.getNumSolutions()):
2065                           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
2066                          # flux=C-util.reorderComponents(B,[1,0])
2067                     else:                     else:
2068                        for l in range(self.getDim()): Z2+=(C[l]-B[l])**2                        for l in range(self.getDim()): flux2+=(C[l]-B[l])**2
2069                  length_of_Z=util.sqrt(Z2)                        #flux=C-B
2070                    length_of_flux=util.sqrt(flux2)
2071              elif C.isEmpty():              elif C.isEmpty():
2072                length_of_Z=util.length(B)                length_of_flux=util.length(B)
2073                  #flux=B
2074              else:              else:
2075                length_of_Z=util.length(C)                length_of_flux=util.length(C)
2076                  #flux=C
2077    
2078              Z_max=util.Lsup(length_of_Z)              #length_of_flux=util.length(flux)
2079              if Z_max>0.:              flux_max=util.Lsup(length_of_flux)
2080                if flux_max>0.:
2081                   # length_of_A=util.inner(flux,util.tensormutiply(A,flux))
2082                 length_of_A=util.length(A)                 length_of_A=util.length(A)
2083                 A_max=util.Lsup(length_of_A)                 A_max=util.Lsup(length_of_A)
2084                 if A_max>0:                 if A_max>0:
2085                      inv_A=1./(length_of_A+A_max*self.__TOL)                      inv_A=1./(length_of_A+A_max*self.__TOL)
2086                 else:                 else:
2087                      inv_A=1./self.__TOL                      inv_A=1./self.__TOL
2088                 peclet_number=length_of_Z*h/2*inv_A                 peclet_number=length_of_flux*h/2*inv_A
2089                 xi=self.__xi(peclet_number)                 xi=self.__xi(peclet_number)
2090                 self.__Xi=h*xi/(length_of_Z+Z_max*self.__TOL)                 self.__Xi=h*xi/(length_of_flux+flux_max*self.__TOL)
2091                 self.trace("preclet number = %e"%util.Lsup(peclet_number))                 self.trace("preclet number = %e"%util.Lsup(peclet_number))
2092        return self.__Xi        return self.__Xi
2093    
# Line 2103  class AdvectivePDE(LinearPDE): Line 2125  class AdvectivePDE(LinearPDE):
2125                        for k in range(self.getNumSolutions()):                        for k in range(self.getNumSolutions()):
2126                           for l in range(self.getDim()):                           for l in range(self.getDim()):
2127                              if not C.isEmpty() and not B.isEmpty():                              if not C.isEmpty() and not B.isEmpty():
2128                                   # tmp=C-util.reorderComponents(B,[0,2,1])
2129                                   # Aout=Aout+Xi*util.generalTensorProduct(util.reorder(tmp,[1,2,0]),tmp,offset=1)
2130                                 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])
2131                              elif C.isEmpty():                              elif C.isEmpty():
2132                                 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]
2133                                   # Aout=Aout+Xi*util.generalTensorProduct(util.reorder(B,[2,1,0]),util.reorder(B,[0,2,1]),offset=1)
2134                              else:                              else:
2135                                 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]
2136                                   # Aout=Aout+Xi*util.generalTensorProduct(util.reorder(C,[1,2,0]),C,offset=1)
2137              else:              else:
2138                  for j in range(self.getDim()):                  for j in range(self.getDim()):
2139                     for l in range(self.getDim()):                     for l in range(self.getDim()):
# Line 2117  class AdvectivePDE(LinearPDE): Line 2143  class AdvectivePDE(LinearPDE):
2143                            Aout[j,l]+=Xi*B[j]*B[l]                            Aout[j,l]+=Xi*B[j]*B[l]
2144                        else:                        else:
2145                            Aout[j,l]+=Xi*C[j]*C[l]                            Aout[j,l]+=Xi*C[j]*C[l]
2146                     # if not C.isEmpty() and not B.isEmpty():
2147                     #    tmp=C-B
2148                     #    Aout=Aout+Xi*util.outer(tmp,tmp)
2149                     # elif C.isEmpty():
2150                     #    Aout=Aout+Xi*util.outer(B,B)
2151                     # else:
2152                     # Aout=Aout+Xi*util.outer(C,C)
2153           return Aout           return Aout
2154       elif name == "B" :       elif name == "B" :
2155           B=self.getCoefficient("B")           B=self.getCoefficient("B")
# Line 2137  class AdvectivePDE(LinearPDE): Line 2170  class AdvectivePDE(LinearPDE):
2170                       for i in range(self.getNumEquations()):                       for i in range(self.getNumEquations()):
2171                          for j in range(self.getDim()):                          for j in range(self.getDim()):
2172                             Bout[i,j,k]+=tmp*C[p,i,j]                             Bout[i,j,k]+=tmp*C[p,i,j]
2173                               # Bout=Bout+Xi*util.generalTensorProduct(util.reorder(C,[1,2,0]),D,offset=1)
2174              else:              else:
2175                 tmp=Xi*D                 tmp=Xi*D
2176                 for j in range(self.getDim()): Bout[j]+=tmp*C[j]                 for j in range(self.getDim()): Bout[j]+=tmp*C[j]
2177                   # Bout=Bout+Xi*D*C
2178           return Bout           return Bout
2179       elif name == "C" :       elif name == "C" :
2180           B=self.getCoefficient("B")           B=self.getCoefficient("B")
# Line 2160  class AdvectivePDE(LinearPDE): Line 2195  class AdvectivePDE(LinearPDE):
2195                        for i in range(self.getNumEquations()):                        for i in range(self.getNumEquations()):
2196                          for l in range(self.getDim()):                          for l in range(self.getDim()):
2197                                   Cout[i,k,l]+=tmp*B[p,l,i]                                   Cout[i,k,l]+=tmp*B[p,l,i]
2198                                     # Cout=Cout+Xi*B[p,l,i]*D[p,k]
2199              else:              else:
2200                 tmp=Xi*D                 tmp=Xi*D
2201                 for j in range(self.getDim()): Cout[j]+=tmp*B[j]                 for j in range(self.getDim()): Cout[j]+=tmp*B[j]
2202                   # Cout=Cout+tmp*D*B
2203           return Cout           return Cout
2204       elif name == "D" :       elif name == "D" :
2205           return self.getCoefficient("D")           return self.getCoefficient("D")
# Line 2186  class AdvectivePDE(LinearPDE): Line 2223  class AdvectivePDE(LinearPDE):
2223                         for j in range(self.getDim()):                         for j in range(self.getDim()):
2224                            if not C.isEmpty() and not B.isEmpty():                            if not C.isEmpty() and not B.isEmpty():
2225                               Xout[i,j]+=tmp*(C[p,i,j]-B[p,j,i])                               Xout[i,j]+=tmp*(C[p,i,j]-B[p,j,i])
2226                                 # Xout=X_out+Xi*util.inner(Y,C-util.reorderComponents(B,[0,2,1]),offset=1)
2227                            elif C.isEmpty():                            elif C.isEmpty():
2228                               Xout[i,j]-=tmp*B[p,j,i]                               Xout[i,j]-=tmp*B[p,j,i]
2229                                 # Xout=X_out-Xi*util.inner(Y,util.reorderComponents(B,[0,2,1]),offset=1)
2230                            else:                            else:
2231                               Xout[i,j]+=tmp*C[p,i,j]                               Xout[i,j]+=tmp*C[p,i,j]
2232                                 # Xout=X_out+Xi*util.inner(Y,C,offset=1)
2233              else:              else:
2234                   tmp=Xi*Y                   tmp=Xi*Y
2235                   for j in range(self.getDim()):                   for j in range(self.getDim()):
2236                      if not C.isEmpty() and not B.isEmpty():                      if not C.isEmpty() and not B.isEmpty():
2237                         Xout[j]+=tmp*(C[j]-B[j])                         Xout[j]+=tmp*(C[j]-B[j])
2238                           # Xout=Xout+Xi*Y*(C-B)
2239                      elif C.isEmpty():                      elif C.isEmpty():
2240                         Xout[j]-=tmp*B[j]                         Xout[j]-=tmp*B[j]
2241                           # Xout=Xout-Xi*Y*B
2242                      else:                      else:
2243                         Xout[j]+=tmp*C[j]                         Xout[j]+=tmp*C[j]
2244                           # Xout=Xout+Xi*Y*C
2245           return Xout           return Xout
2246       elif name == "Y" :       elif name == "Y" :
2247           return self.getCoefficient("Y")           return self.getCoefficient("Y")
# Line 2213  class AdvectivePDE(LinearPDE): Line 2256  class AdvectivePDE(LinearPDE):
2256       elif name == "r" :       elif name == "r" :
2257           return self.getCoefficient("r")           return self.getCoefficient("r")
2258       elif name == "q" :       elif name == "q" :
2259             return self.getCoefficient("q")
2260         else:
2261            raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
2262    
2263    class AdvectionDiffusion(LinearPDE):
2264       """
2265       Class to define PDE equation of the unisotropic advection-diffusion problem, which is genear L{LinearPDE} of the form
2266    
2267       M{S{omega}*u + inner(v,grad(u))- grad(matrixmult(k_bar,grad(u))[j])[j] = f}
2268    
2269       with natural boundary conditons
2270    
2271       M{n[j]*matrixmult(k,grad(u))[j] = g- S{alpha}u }
2272    
2273       and constraints:
2274    
2275       M{u=r} where M{q>0}
2276    
2277       and
2278    
2279       M{k_bar[i,j]=k[i,j]+upwind[i]*upwind[j]}
2280    
2281       """
2282    
2283       def __init__(self,domain,debug=False):
2284         """
2285         initializes a new Poisson equation
2286    
2287         @param domain: domain of the PDE
2288         @type domain: L{Domain<escript.Domain>}
2289         @param debug: if True debug informations are printed.
2290    
2291         """
2292         super(AdvectionDiffusion, self).__init__(domain,1,1,debug)
2293         self.COEFFICIENTS={"omega": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
2294                            "k": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_DIM,PDECoefficient.BY_DIM),PDECoefficient.OPERATOR),
2295                            "f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
2296                            "v": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_DIM,),PDECoefficient.OPERATOR),
2297                            "upwind": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_DIM,),PDECoefficient.OPERATOR),
2298                            "alpha": PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
2299                            "g": PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
2300                            "r": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH),
2301                            "q": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}
2302    
2303       def setValue(self,**coefficients):
2304         """
2305         sets new values to coefficients
2306    
2307         @param coefficients: new values assigned to coefficients
2308         @keyword omega: value for coefficient M{S{omega}}
2309         @type omega: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
2310         @keyword k: value for coefficient M{k}
2311         @type k: any type that can be casted to L{Tensor<escript.Tensor>} object on L{Function<escript.Function>}.
2312         @keyword v: value for coefficient M{v}
2313         @type v: any type that can be casted to L{Vector<escript.Vector>} object on L{Function<escript.Function>}.
2314         @keyword upwind: value for upwind term M{upwind}
2315         @type upwind: any type that can be casted to L{Vector<escript.Vector>} object on L{Function<escript.Function>}.
2316         @keyword f: value for right hand side M{f}
2317         @type f: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
2318         @keyword alpha: value for right hand side M{S{alpha}}
2319         @type alpha: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
2320         @keyword g: value for right hand side M{g}
2321         @type g: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
2322         @keyword r: prescribed values M{r} for the solution in constraints.
2323         @type r: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2324                   depending of reduced order is used for the representation of the equation.
2325         @keyword q: mask for location of constraints
2326         @type q: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2327                   depending of reduced order is used for the representation of the equation.
2328         @raise IllegalCoefficient: if an unknown coefficient keyword is used.
2329         """
2330         super(AdvectionDiffusion, self).setValue(**coefficients)
2331    
2332       def getCoefficientOfGeneralPDE(self,name):
2333         """
2334         return the value of the coefficient name of the general PDE
2335    
2336         @param name: name of the coefficient requested.
2337         @type name: C{string}
2338         @return: the value of the coefficient  name
2339         @rtype: L{Data<escript.Data>}
2340         @raise IllegalCoefficient: if name is not one of coefficients
2341                      "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}.
2342         @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
2343         """
2344         if name == "A" :
2345             return self.getCoefficient("k")+util.outer(self.getCoefficient("upwind"),self.getCoefficient("upwind"))
2346         elif name == "B" :
2347             return escript.Data()
2348         elif name == "C" :
2349             return self.getCoefficient("v")
2350         elif name == "D" :
2351             return self.getCoefficient("omega")
2352         elif name == "X" :
2353             return escript.Data()
2354         elif name == "Y" :
2355             return self.getCoefficient("f")
2356         elif name == "d" :
2357             return self.getCoefficient("alpha")
2358         elif name == "y" :
2359             return self.getCoefficient("g")
2360         elif name == "d_contact" :
2361             return escript.Data()
2362         elif name == "y_contact" :
2363             return escript.Data()
2364         elif name == "r" :
2365             return self.getCoefficient("r")
2366         elif name == "q" :
2367           return self.getCoefficient("q")           return self.getCoefficient("q")
2368       else:       else:
2369          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.430

  ViewVC Help
Powered by ViewVC 1.1.26