/[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 410 by gross, Fri Dec 23 01:27:38 2005 UTC revision 969 by ksteube, Tue Feb 13 23:02:23 2007 UTC
# Line 1  Line 1 
1  # $Id$  # $Id$
   
 #  
 #      COPYRIGHT ACcESS 2004 -  All Rights Reserved  
 #  
 #   This software is the property of ACcESS.  No part of this code  
 #   may be copied in any form or by any means without the expressed written  
 #   consent of ACcESS.  Copying, use or modification of this software  
 #   by any unauthorised person is illegal unless that  
 #   person has a software license agreement with ACcESS.  
 #  
2  """  """
3  The module provides an interface to define and solve linear partial  The module provides an interface to define and solve linear partial
4  differential equations (PDEs) within L{escript}. L{linearPDEs} does not provide any  differential equations (PDEs) within L{escript}. L{linearPDEs} does not provide any
# Line 17  the PDE solver library defined through t Line 7  the PDE solver library defined through t
7  The general interface is provided through the L{LinearPDE} class. The  The general interface is provided through the L{LinearPDE} class. The
8  L{AdvectivePDE} which is derived from the L{LinearPDE} class  L{AdvectivePDE} which is derived from the L{LinearPDE} class
9  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},
10  L{Helmholtz}, L{LameEquation}, L{AdvectionDiffusion}  L{Helmholtz}, L{LameEquation}, L{AdvectivePDE}
11  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
12  to define of solve these sepecial PDEs.  to define of solve these sepecial PDEs.
13    
14  @var __author__: name of author  @var __author__: name of author
15  @var __licence__: licence agreement  @var __copyright__: copyrights
16    @var __license__: licence agreement
17  @var __url__: url entry point on documentation  @var __url__: url entry point on documentation
18  @var __version__: version  @var __version__: version
19  @var __date__: date of the version  @var __date__: date of the version
# Line 33  import util Line 24  import util
24  import numarray  import numarray
25    
26  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
27  __licence__="contact: esys@access.uq.edu.au"  __copyright__="""  Copyright (c) 2006 by ACcESS MNRF
28  __url__="http://www.iservo.edu.au/esys/escript"                      http://www.access.edu.au
29                    Primary Business: Queensland, Australia"""
30    __license__="""Licensed under the Open Software License version 3.0
31                 http://www.opensource.org/licenses/osl-3.0.php"""
32    __url__="http://www.iservo.edu.au/esys"
33  __version__="$Revision$"  __version__="$Revision$"
34  __date__="$Date$"  __date__="$Date$"
35    
# Line 120  class PDECoefficient(object): Line 115  class PDECoefficient(object):
115         @param domain: domain on which the PDE uses the coefficient         @param domain: domain on which the PDE uses the coefficient
116         @type domain: L{Domain<escript.Domain>}         @type domain: L{Domain<escript.Domain>}
117         @param reducedEquationOrder: True to indicate that reduced order is used to represent the equation         @param reducedEquationOrder: True to indicate that reduced order is used to represent the equation
118         @type domain: C{bool}         @type reducedEquationOrder: C{bool}
119         @param reducedSolutionOrder: True to indicate that reduced order is used to represent the solution         @param reducedSolutionOrder: True to indicate that reduced order is used to represent the solution
120         @type domain: C{bool}         @type reducedSolutionOrder: C{bool}
121         @return:  L{FunctionSpace<escript.FunctionSpace>} of the coefficient         @return:  L{FunctionSpace<escript.FunctionSpace>} of the coefficient
122         @rtype:  L{FunctionSpace<escript.FunctionSpace>}         @rtype:  L{FunctionSpace<escript.FunctionSpace>}
123         """         """
# Line 160  class PDECoefficient(object): Line 155  class PDECoefficient(object):
155         @param numSolutions: number of components of the PDE solution         @param numSolutions: number of components of the PDE solution
156         @type numSolutions: C{int}         @type numSolutions: C{int}
157         @param reducedEquationOrder: True to indicate that reduced order is used to represent the equation         @param reducedEquationOrder: True to indicate that reduced order is used to represent the equation
158         @type domain: C{bool}         @type reducedEquationOrder: C{bool}
159         @param reducedSolutionOrder: True to indicate that reduced order is used to represent the solution         @param reducedSolutionOrder: True to indicate that reduced order is used to represent the solution
160         @type domain: C{bool}         @type reducedSolutionOrder: C{bool}
161         @param newValue: number of components of the PDE solution         @param newValue: number of components of the PDE solution
162         @type newValue: any object that can be converted into a L{Data<escript.Data>} object with the appropriate shape and L{FunctionSpace<escript.FunctionSpace>}         @type newValue: any object that can be converted into a L{Data<escript.Data>} object with the appropriate shape and L{FunctionSpace<escript.FunctionSpace>}
163         @raise IllegalCoefficientValue: if the shape of the assigned value does not match the shape of the coefficient         @raise IllegalCoefficientValue: if the shape of the assigned value does not match the shape of the coefficient
# Line 418  class LinearPDE(object): Line 413  class LinearPDE(object):
413     @cvar SCSL: SGI SCSL solver library     @cvar SCSL: SGI SCSL solver library
414     @cvar MKL: Intel's MKL solver library     @cvar MKL: Intel's MKL solver library
415     @cvar UMFPACK: the UMFPACK library     @cvar UMFPACK: the UMFPACK library
416       @cvar TRILINOS: the TRILINOS parallel solver class library from Sandia Natl Labs
417     @cvar ITERATIVE: The default iterative solver     @cvar ITERATIVE: The default iterative solver
418       @cvar AMG: algebraic multi grid
419       @cvar RILU: recursive ILU
420    
421     """     """
422     DEFAULT= 0     DEFAULT= 0
# Line 443  class LinearPDE(object): Line 441  class LinearPDE(object):
441     UMFPACK= 16     UMFPACK= 16
442     ITERATIVE= 20     ITERATIVE= 20
443     PASO= 21     PASO= 21
444       AMG= 22
445       RILU = 23
446       TRILINOS = 24
447    
448     __TOL=1.e-13     SMALL_TOLERANCE=1.e-13
449     __PACKAGE_KEY="package"     __PACKAGE_KEY="package"
450     __METHOD_KEY="method"     __METHOD_KEY="method"
451     __SYMMETRY_KEY="symmetric"     __SYMMETRY_KEY="symmetric"
# Line 665  class LinearPDE(object): Line 666  class LinearPDE(object):
666       @rtype: L{Data<escript.Data>}       @rtype: L{Data<escript.Data>}
667       """       """
668       if u==None:       if u==None:
669            return self.getOperator()*self.getSolution()          return self.getOperator()*self.getSolution()
670       else:       else:
671          self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())          return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())
672    
673     def getResidual(self,u=None):     def getResidual(self,u=None):
674       """       """
# Line 699  class LinearPDE(object): Line 700  class LinearPDE(object):
700        else:        else:
701           A=self.getCoefficientOfGeneralPDE("A")           A=self.getCoefficientOfGeneralPDE("A")
702           if not A.isEmpty():           if not A.isEmpty():
703              tol=util.Lsup(A)*self.__TOL              tol=util.Lsup(A)*self.SMALL_TOLERANCE
704              if self.getNumSolutions()>1:              if self.getNumSolutions()>1:
705                 for i in range(self.getNumEquations()):                 for i in range(self.getNumEquations()):
706                    for j in range(self.getDim()):                    for j in range(self.getDim()):
# Line 723  class LinearPDE(object): Line 724  class LinearPDE(object):
724              if verbose: print "non-symmetric PDE because C is not present but B is"              if verbose: print "non-symmetric PDE because C is not present but B is"
725              out=False              out=False
726           elif not B.isEmpty() and not C.isEmpty():           elif not B.isEmpty() and not C.isEmpty():
727              tol=(util.Lsup(B)+util.Lsup(C))*self.__TOL/2.              tol=(util.Lsup(B)+util.Lsup(C))*self.SMALL_TOLERANCE/2.
728              if self.getNumSolutions()>1:              if self.getNumSolutions()>1:
729                 for i in range(self.getNumEquations()):                 for i in range(self.getNumEquations()):
730                     for j in range(self.getDim()):                     for j in range(self.getDim()):
# Line 739  class LinearPDE(object): Line 740  class LinearPDE(object):
740           if self.getNumSolutions()>1:           if self.getNumSolutions()>1:
741             D=self.getCoefficientOfGeneralPDE("D")             D=self.getCoefficientOfGeneralPDE("D")
742             if not D.isEmpty():             if not D.isEmpty():
743               tol=util.Lsup(D)*self.__TOL               tol=util.Lsup(D)*self.SMALL_TOLERANCE
744               for i in range(self.getNumEquations()):               for i in range(self.getNumEquations()):
745                  for k in range(self.getNumSolutions()):                  for k in range(self.getNumSolutions()):
746                    if util.Lsup(D[i,k]-D[k,i])>tol:                    if util.Lsup(D[i,k]-D[k,i])>tol:
# Line 747  class LinearPDE(object): Line 748  class LinearPDE(object):
748                        out=False                        out=False
749             d=self.getCoefficientOfGeneralPDE("d")             d=self.getCoefficientOfGeneralPDE("d")
750             if not d.isEmpty():             if not d.isEmpty():
751               tol=util.Lsup(d)*self.__TOL               tol=util.Lsup(d)*self.SMALL_TOLERANCE
752               for i in range(self.getNumEquations()):               for i in range(self.getNumEquations()):
753                  for k in range(self.getNumSolutions()):                  for k in range(self.getNumSolutions()):
754                    if util.Lsup(d[i,k]-d[k,i])>tol:                    if util.Lsup(d[i,k]-d[k,i])>tol:
# Line 755  class LinearPDE(object): Line 756  class LinearPDE(object):
756                        out=False                        out=False
757             d_contact=self.getCoefficientOfGeneralPDE("d_contact")             d_contact=self.getCoefficientOfGeneralPDE("d_contact")
758             if not d_contact.isEmpty():             if not d_contact.isEmpty():
759               tol=util.Lsup(d_contact)*self.__TOL               tol=util.Lsup(d_contact)*self.SMALL_TOLERANCE
760               for i in range(self.getNumEquations()):               for i in range(self.getNumEquations()):
761                  for k in range(self.getNumSolutions()):                  for k in range(self.getNumSolutions()):
762                    if util.Lsup(d_contact[i,k]-d_contact[k,i])>tol:                    if util.Lsup(d_contact[i,k]-d_contact[k,i])>tol:
# Line 821  class LinearPDE(object): Line 822  class LinearPDE(object):
822         sets a new solver         sets a new solver
823    
824         @param solver: sets a new solver method.         @param solver: sets a new solver method.
825         @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}, L{AMG}
826         @param preconditioner: sets a new solver method.         @param preconditioner: sets a new solver method.
827         @type solver: one of L{DEFAULT}, L{JACOBI} L{ILU0}, L{ILUT},L{SSOR}         @type preconditioner: one of L{DEFAULT}, L{JACOBI} L{ILU0}, L{ILUT},L{SSOR}, L{RILU}
828         """         """
829         if solver==None: solve=self.DEFAULT         if solver==None: solve=self.DEFAULT
830         if preconditioner==None: preconditioner=self.DEFAULT         if preconditioner==None: preconditioner=self.DEFAULT
# Line 843  class LinearPDE(object): Line 844  class LinearPDE(object):
844    
845         m=self.getSolverMethod()         m=self.getSolverMethod()
846         p=self.getSolverPackage()         p=self.getSolverPackage()
847           method=""
848         if m[0]==self.DEFAULT: method="DEFAULT"         if m[0]==self.DEFAULT: method="DEFAULT"
849         elif m[0]==self.DIRECT: method= "DIRECT"         elif m[0]==self.DIRECT: method= "DIRECT"
850         elif m[0]==self.ITERATIVE: method= "ITERATIVE"         elif m[0]==self.ITERATIVE: method= "ITERATIVE"
# Line 855  class LinearPDE(object): Line 857  class LinearPDE(object):
857         elif m[0]==self.GMRES: method= "GMRES"         elif m[0]==self.GMRES: method= "GMRES"
858         elif m[0]==self.PRES20: method= "PRES20"         elif m[0]==self.PRES20: method= "PRES20"
859         elif m[0]==self.LUMPING: method= "LUMPING"         elif m[0]==self.LUMPING: method= "LUMPING"
860         else : method="unknown"         elif m[0]==self.AMG: method= "AMG"
861         if m[1]==self.DEFAULT: method+="DEFAULT"         if m[1]==self.DEFAULT: method+="+DEFAULT"
862         elif m[1]==self.JACOBI: method+= "JACOBI"         elif m[1]==self.JACOBI: method+= "+JACOBI"
863         elif m[1]==self.ILU0: method+= "ILU0"         elif m[1]==self.ILU0: method+= "+ILU0"
864         elif m[1]==self.ILUT: method+= "ILUT"         elif m[1]==self.ILUT: method+= "+ILUT"
865         elif m[1]==self.SSOR: method+= "SSOR"         elif m[1]==self.SSOR: method+= "+SSOR"
866         else : method+="unknown"         elif m[1]==self.AMG: method+= "+AMG"
867           elif m[1]==self.RILU: method+= "+RILU"
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"
871         elif p==self.SCSL: package= "SCSL"         elif p==self.SCSL: package= "SCSL"
872         elif p==self.UMFPACK: package= "UMFPACK"         elif p==self.UMFPACK: package= "UMFPACK"
873           elif p==self.TRILINOS: package= "TRILINOS"
874         else : method="unknown"         else : method="unknown"
875         return "%s solver of %s package"%(method,package)         return "%s solver of %s package"%(method,package)
876    
# Line 884  class LinearPDE(object): Line 888  class LinearPDE(object):
888         """         """
889         sets a new solver package         sets a new solver package
890    
891         @param solver: sets a new solver method.         @param package: sets a new solver method.
892         @type solver: one of L{DEFAULT}, L{PASO} L{SCSL}, L{MKL}, L{UMLPACK}         @type package: one of L{DEFAULT}, L{PASO} L{SCSL}, L{MKL}, L{UMFPACK}, L{TRILINOS}
893         """         """
894         if package==None: package=self.DEFAULT         if package==None: package=self.DEFAULT
895         if not package==self.getSolverPackage():         if not package==self.getSolverPackage():
896             self.__solver_method=solver             self.__solver_package=package
897             self.__checkMatrixType()             self.__checkMatrixType()
898             self.trace("New solver is %s"%self.getSolverMethodName())             self.trace("New solver is %s"%self.getSolverMethodName())
899    
# Line 922  class LinearPDE(object): Line 926  class LinearPDE(object):
926         @param tol: new tolerance for the solver. If the tol is lower then the current tolerence         @param tol: new tolerance for the solver. If the tol is lower then the current tolerence
927                     the system will be resolved.                     the system will be resolved.
928         @type tol: positive C{float}         @type tol: positive C{float}
929         @raise ValueException: if tolerance is not positive.         @raise ValueError: if tolerance is not positive.
930         """         """
931         if not tol>0:         if not tol>0:
932             raise ValueException,"Tolerance as to be positive"             raise ValueError,"Tolerance as to be positive"
933         if tol<self.getTolerance(): self.__invalidateSolution()         if tol<self.getTolerance(): self.__invalidateSolution()
934         self.trace("New tolerance %e"%tol)         self.trace("New tolerance %e"%tol)
935         self.__tolerance=tol         self.__tolerance=tol
# Line 1512  class LinearPDE(object): Line 1516  class LinearPDE(object):
1516           r=self.getCoefficientOfGeneralPDE("r")           r=self.getCoefficientOfGeneralPDE("r")
1517           homogeneous_constraint=True           homogeneous_constraint=True
1518           if not q.isEmpty() and not r.isEmpty():           if not q.isEmpty() and not r.isEmpty():
1519               if util.Lsup(q*r)>=1.e-13*util.Lsup(r):               if util.Lsup(q*r)>0.:
1520                 self.trace("Inhomogeneous constraint detected.")                 self.trace("Inhomogeneous constraint detected.")
1521                 self.__invalidateSystem()                 self.__invalidateSystem()
1522    
# Line 1526  class LinearPDE(object): Line 1530  class LinearPDE(object):
1530         if not self.__operator_is_Valid or not self.__righthandside_isValid:         if not self.__operator_is_Valid or not self.__righthandside_isValid:
1531            if self.isUsingLumping():            if self.isUsingLumping():
1532                if not self.__operator_is_Valid:                if not self.__operator_is_Valid:
1533                   if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution(): raise TypeError,"Lumped matrix requires same order for equations and unknowns"                   if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():
1534                   if not self.getCoefficientOfGeneralPDE("A").isEmpty(): raise Warning,"Using coefficient A in lumped matrix can produce wrong results"                        raise TypeError,"Lumped matrix requires same order for equations and unknowns"
1535                   if not self.getCoefficientOfGeneralPDE("B").isEmpty(): raise Warning,"Using coefficient B in lumped matrix can produce wrong results"                   if not self.getCoefficientOfGeneralPDE("A").isEmpty():
1536                   if not self.getCoefficientOfGeneralPDE("C").isEmpty(): raise Warning,"Using coefficient C in lumped matrix can produce wrong results"                        raise ValueError,"coefficient A in lumped matrix may not be present."
1537                   mat=self.__getNewOperator()                   if not self.getCoefficientOfGeneralPDE("B").isEmpty():
1538                   self.getDomain().addPDEToSystem(mat,escript.Data(), \                        raise ValueError,"coefficient A in lumped matrix may not be present."
1539                             self.getCoefficientOfGeneralPDE("A"), \                   if not self.getCoefficientOfGeneralPDE("C").isEmpty():
1540                             self.getCoefficientOfGeneralPDE("B"), \                        raise ValueError,"coefficient A in lumped matrix may not be present."
1541                             self.getCoefficientOfGeneralPDE("C"), \                   D=self.getCoefficientOfGeneralPDE("D")
1542                             self.getCoefficientOfGeneralPDE("D"), \                   if not D.isEmpty():
1543                             escript.Data(), \                       if self.getNumSolutions()>1:
1544                             escript.Data(), \                          #D_times_e=util.matrixmult(D,numarray.ones((self.getNumSolutions(),)))
1545                             self.getCoefficientOfGeneralPDE("d"), \                          D_times_e=util.matrix_mult(D,numarray.ones((self.getNumSolutions(),)))
1546                             escript.Data(),\                       else:
1547                             self.getCoefficientOfGeneralPDE("d_contact"), \                          D_times_e=D
1548                             escript.Data())                   else:
1549                   self.__operator=1./(mat*escript.Data(1,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True))                      D_times_e=escript.Data()
1550                   del mat                   d=self.getCoefficientOfGeneralPDE("d")
1551                     if not d.isEmpty():
1552                         if self.getNumSolutions()>1:
1553                            #d_times_e=util.matrixmult(d,numarray.ones((self.getNumSolutions(),)))
1554                            d_times_e=util.matrix_mult(d,numarray.ones((self.getNumSolutions(),)))
1555                         else:
1556                            d_times_e=d
1557                     else:
1558                        d_times_e=escript.Data()
1559                     d_contact=self.getCoefficientOfGeneralPDE("d_contact")
1560                     if not d_contact.isEmpty():
1561                         if self.getNumSolutions()>1:
1562                            d_contact_times_e=util.matrixmult(d_contact,numarray.ones((self.getNumSolutions(),)))
1563                         else:
1564                            d_contact_times_e=d_contact
1565                     else:
1566                        d_contact_times_e=escript.Data()
1567        
1568                     self.__operator=self.__getNewRightHandSide()
1569                     self.getDomain().addPDEToRHS(self.__operator, \
1570                                                  escript.Data(), \
1571                                                  D_times_e, \
1572                                                  d_times_e,\
1573                                                  d_contact_times_e)
1574                     self.__operator=1./self.__operator
1575                   self.trace("New lumped operator has been built.")                   self.trace("New lumped operator has been built.")
1576                   self.__operator_is_Valid=True                   self.__operator_is_Valid=True
1577                if not self.__righthandside_isValid:                if not self.__righthandside_isValid:
# Line 1784  class LameEquation(LinearPDE): Line 1812  class LameEquation(LinearPDE):
1812     """     """
1813     Class to define a Lame equation problem:     Class to define a Lame equation problem:
1814    
1815     M{-grad(S{mu}*(grad(u[i])[j]+grad(u[j])[i]))[j] - grad(S{lambda}*grad(u[j])[i])[j] = F_i -grad(S{sigma}[i,j])[j] }     M{-grad(S{mu}*(grad(u[i])[j]+grad(u[j])[i]))[j] - grad(S{lambda}*grad(u[k])[k])[j] = F_i -grad(S{sigma}[ij])[j] }
1816    
1817     with natural boundary conditons:     with natural boundary conditons:
1818    
1819     M{n[j]*(S{mu}*(grad(u[i])[j]+grad(u[j])[i]) - S{lambda}*grad(u[j])[i]) = f_i -n[j]*S{sigma}[i,j] }     M{n[j]*(S{mu}*(grad(u[i])[j]+grad(u[j])[i]) + S{lambda}*grad(u[k])[k]) = f_i +n[j]*S{sigma}[ij] }
1820    
1821     and constraints:     and constraints:
1822    
# Line 1808  class LameEquation(LinearPDE): Line 1836  class LameEquation(LinearPDE):
1836                            "q"            : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}                            "q"            : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}
1837        self.setSymmetryOn()        self.setSymmetryOn()
1838    
1839     def setValue(self,**coefficients):     def setValues(self,**coefficients):
1840       """       """
1841       sets new values to coefficients       sets new values to coefficients
1842    
# Line 1831  class LameEquation(LinearPDE): Line 1859  class LameEquation(LinearPDE):
1859                 depending of reduced order is used for the representation of the equation.                 depending of reduced order is used for the representation of the equation.
1860       @raise IllegalCoefficient: if an unknown coefficient keyword is used.       @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1861       """       """
1862       super(LameEquation, self).setValue(**coefficients)       super(LameEquation, self).setValues(**coefficients)
1863    
1864     def getCoefficientOfGeneralPDE(self,name):     def getCoefficientOfGeneralPDE(self,name):
1865       """       """
# Line 1945  class AdvectivePDE(LinearPDE): Line 1973  class AdvectivePDE(LinearPDE):
1973           self.__xi=xi           self.__xi=xi
1974        self.__Xi=escript.Data()        self.__Xi=escript.Data()
1975    
1976     def setValue(**coefficients):     def setValue(self,**coefficients):
1977        """        """
1978        sets new values to coefficients        sets new values to coefficients
1979    
# Line 2029  class AdvectivePDE(LinearPDE): Line 2057  class AdvectivePDE(LinearPDE):
2057       """       """
2058       return escript.Scalar(0.5,P.getFunctionSpace())       return escript.Scalar(0.5,P.getFunctionSpace())
2059    
    def __calculateXi(self,peclet_factor,flux,h):  
        flux=util.Lsup(flux)  
        if flux_max>0.:  
           return h*self.__xi(flux*peclet_factor)/(flux+flux_max*self.__TOL)  
        else:  
           return 0.  
   
2060     def __getXi(self):     def __getXi(self):
2061        if self.__Xi.isEmpty():        if self.__Xi.isEmpty():
2062           B=self.getCoefficient("B")           B=self.getCoefficient("B")
# Line 2045  class AdvectivePDE(LinearPDE): Line 2066  class AdvectivePDE(LinearPDE):
2066           self.__Xi=escript.Scalar(0.,self.getFunctionSpaceForCoefficient("A"))           self.__Xi=escript.Scalar(0.,self.getFunctionSpaceForCoefficient("A"))
2067           if not C.isEmpty() or not B.isEmpty():           if not C.isEmpty() or not B.isEmpty():
2068              if not C.isEmpty() and not B.isEmpty():              if not C.isEmpty() and not B.isEmpty():
                 flux2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))  
2069                  if self.getNumEquations()>1:                  if self.getNumEquations()>1:
2070                     if self.getNumSolutions()>1:                     if self.getNumSolutions()>1:
2071                          flux2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))
2072                        for i in range(self.getNumEquations()):                        for i in range(self.getNumEquations()):
2073                           for k in range(self.getNumSolutions()):                           for k in range(self.getNumSolutions()):
2074                              for l in range(self.getDim()): flux2+=(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
2075                          length_of_flux=util.sqrt(flux2)
2076                        # flux=C-util.reorderComponents(B,[0,2,1])                        # flux=C-util.reorderComponents(B,[0,2,1])
2077                     else:                     else:
2078                          flux2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))
2079                        for i in range(self.getNumEquations()):                        for i in range(self.getNumEquations()):
2080                           for l in range(self.getDim()): flux2+=(C[i,l]-B[i,l])**2                           for l in range(self.getDim()): flux2+=(C[i,l]-B[i,l])**2
2081                          length_of_flux=util.sqrt(flux2)
2082                        # flux=C-B                        # flux=C-B
2083                  else:                  else:
2084                     if self.getNumSolutions()>1:                     if self.getNumSolutions()>1:
2085                          flux2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))
2086                        for k in range(self.getNumSolutions()):                        for k in range(self.getNumSolutions()):
2087                           for l in range(self.getDim()): flux2+=(C[k,l]-B[l,k])**2                           for l in range(self.getDim()): flux2+=(C[k,l]-B[l,k])**2
2088                        # flux=C-util.reorderComponents(B,[1,0])                        # flux=C-util.reorderComponents(B,[1,0])
2089                          length_of_flux=util.sqrt(flux2)
2090                     else:                     else:
2091                        for l in range(self.getDim()): flux2+=(C[l]-B[l])**2                        length_of_flux=util.length(C-B)
                       #flux=C-B  
                 length_of_flux=util.sqrt(flux2)  
2092              elif C.isEmpty():              elif C.isEmpty():
2093                length_of_flux=util.length(B)                length_of_flux=util.length(B)
               #flux=B  
2094              else:              else:
2095                length_of_flux=util.length(C)                length_of_flux=util.length(C)
               #flux=C  
   
             #length_of_flux=util.length(flux)  
2096              flux_max=util.Lsup(length_of_flux)              flux_max=util.Lsup(length_of_flux)
2097              if flux_max>0.:              if flux_max>0.:
2098                 # length_of_A=util.inner(flux,util.tensormutiply(A,flux))                if A.isEmpty():
2099                 length_of_A=util.length(A)                    inv_A=1./self.SMALL_TOLERANCE
2100                 A_max=util.Lsup(length_of_A)                    peclet_number=escript.Scalar(inv_A,length_of_flux.getFunctionSpace())
2101                 if A_max>0:                    xi=self.__xi(self,peclet_number)
2102                      inv_A=1./(length_of_A+A_max*self.__TOL)                else:
2103                 else:                    # length_of_A=util.inner(flux,util.tensormutiply(A,flux))
2104                      inv_A=1./self.__TOL                    length_of_A=util.length(A)
2105                 peclet_number=length_of_flux*h/2*inv_A                    A_max=util.Lsup(length_of_A)
2106                 xi=self.__xi(peclet_number)                    if A_max>0:
2107                 self.__Xi=h*xi/(length_of_flux+flux_max*self.__TOL)                         inv_A=1./(length_of_A+A_max*self.SMALL_TOLERANCE)
2108                 self.trace("preclet number = %e"%util.Lsup(peclet_number))                    else:
2109                           inv_A=1./self.SMALL_TOLERANCE
2110                      peclet_number=length_of_flux*h/2*inv_A
2111                      xi=self.__xi(self,peclet_number)
2112                  self.__Xi=h*xi/(length_of_flux+flux_max*self.SMALL_TOLERANCE)
2113                  self.trace("preclet number = %e"%util.Lsup(peclet_number))
2114                else:
2115                  self.__Xi=escript.Scalar(0.,length_of_flux.getFunctionSpace())
2116        return self.__Xi        return self.__Xi
2117    
2118    
# Line 2112  class AdvectivePDE(LinearPDE): Line 2139  class AdvectivePDE(LinearPDE):
2139              Aout=A              Aout=A
2140           else:           else:
2141              if A.isEmpty():              if A.isEmpty():
2142                 Aout=self.createNewCoefficient("A")                 Aout=self.createCoefficientOfGeneralPDE("A")
2143              else:              else:
2144                 Aout=A[:]                 Aout=A[:]
2145              Xi=self.__getXi()              Xi=self.__getXi()
# Line 2132  class AdvectivePDE(LinearPDE): Line 2159  class AdvectivePDE(LinearPDE):
2159                                 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]
2160                                 # Aout=Aout+Xi*util.generalTensorProduct(util.reorder(C,[1,2,0]),C,offset=1)                                 # Aout=Aout+Xi*util.generalTensorProduct(util.reorder(C,[1,2,0]),C,offset=1)
2161              else:              else:
2162                  for j in range(self.getDim()):                 if not C.isEmpty() and not B.isEmpty():
2163                     for l in range(self.getDim()):                     delta=(C-B)
2164                        if not C.isEmpty() and not B.isEmpty():                     Aout+=util.outer(Xi*delta,delta)
2165                            Aout[j,l]+=Xi*(C[j]-B[j])*(C[l]-B[l])                 elif not B.isEmpty():
2166                        elif C.isEmpty():                     Aout+=util.outer(Xi*B,B)
2167                            Aout[j,l]+=Xi*B[j]*B[l]                 elif not C.isEmpty():
2168                        else:                     Aout+=util.outer(Xi*C,C)
                           Aout[j,l]+=Xi*C[j]*C[l]  
                  # if not C.isEmpty() and not B.isEmpty():  
                  #    tmp=C-B  
                  #    Aout=Aout+Xi*util.outer(tmp,tmp)  
                  # elif C.isEmpty():  
                  #    Aout=Aout+Xi*util.outer(B,B)  
                  # else:  
                  # Aout=Aout+Xi*util.outer(C,C)  
2169           return Aout           return Aout
2170       elif name == "B" :       elif name == "B" :
2171             # return self.getCoefficient("B")
2172           B=self.getCoefficient("B")           B=self.getCoefficient("B")
2173           C=self.getCoefficient("C")           C=self.getCoefficient("C")
2174           D=self.getCoefficient("D")           D=self.getCoefficient("D")
# Line 2157  class AdvectivePDE(LinearPDE): Line 2177  class AdvectivePDE(LinearPDE):
2177           else:           else:
2178              Xi=self.__getXi()              Xi=self.__getXi()
2179              if B.isEmpty():              if B.isEmpty():
2180                  Bout=self.createNewCoefficient("B")                  Bout=self.createCoefficientOfGeneralPDE("B")
2181              else:              else:
2182                  Bout=B[:]                  Bout=B[:]
2183              if self.getNumEquations()>1:              if self.getNumEquations()>1:
# Line 2169  class AdvectivePDE(LinearPDE): Line 2189  class AdvectivePDE(LinearPDE):
2189                             Bout[i,j,k]+=tmp*C[p,i,j]                             Bout[i,j,k]+=tmp*C[p,i,j]
2190                             # Bout=Bout+Xi*util.generalTensorProduct(util.reorder(C,[1,2,0]),D,offset=1)                             # Bout=Bout+Xi*util.generalTensorProduct(util.reorder(C,[1,2,0]),D,offset=1)
2191              else:              else:
2192                 tmp=Xi*D                 Bout+=(Xi*D)*C
                for j in range(self.getDim()): Bout[j]+=tmp*C[j]  
                # Bout=Bout+Xi*D*C  
2193           return Bout           return Bout
2194       elif name == "C" :       elif name == "C" :
2195             # return self.getCoefficient("C")
2196           B=self.getCoefficient("B")           B=self.getCoefficient("B")
2197           C=self.getCoefficient("C")           C=self.getCoefficient("C")
2198           D=self.getCoefficient("D")           D=self.getCoefficient("D")
# Line 2182  class AdvectivePDE(LinearPDE): Line 2201  class AdvectivePDE(LinearPDE):
2201           else:           else:
2202              Xi=self.__getXi()              Xi=self.__getXi()
2203              if C.isEmpty():              if C.isEmpty():
2204                  Cout=self.createNewCoefficient("C")                  Cout=self.createCoefficientOfGeneralPDE("C")
2205              else:              else:
2206                  Cout=C[:]                  Cout=C[:]
2207              if self.getNumEquations()>1:              if self.getNumEquations()>1:
# Line 2194  class AdvectivePDE(LinearPDE): Line 2213  class AdvectivePDE(LinearPDE):
2213                                   Cout[i,k,l]+=tmp*B[p,l,i]                                   Cout[i,k,l]+=tmp*B[p,l,i]
2214                                   # Cout=Cout+Xi*B[p,l,i]*D[p,k]                                   # Cout=Cout+Xi*B[p,l,i]*D[p,k]
2215              else:              else:
2216                 tmp=Xi*D                 Cout+=(Xi*D)*B
                for j in range(self.getDim()): Cout[j]+=tmp*B[j]  
                # Cout=Cout+tmp*D*B  
2217           return Cout           return Cout
2218       elif name == "D" :       elif name == "D" :
2219           return self.getCoefficient("D")           return self.getCoefficient("D")
2220       elif name == "X" :       elif name == "X" :
2221             # return self.getCoefficient("X")
2222           X=self.getCoefficient("X")           X=self.getCoefficient("X")
2223           Y=self.getCoefficient("Y")           Y=self.getCoefficient("Y")
2224           B=self.getCoefficient("B")           B=self.getCoefficient("B")
# Line 2209  class AdvectivePDE(LinearPDE): Line 2227  class AdvectivePDE(LinearPDE):
2227              Xout=X              Xout=X
2228           else:           else:
2229              if X.isEmpty():              if X.isEmpty():
2230                  Xout=self.createNewCoefficient("X")                  Xout=self.createCoefficientOfGeneralPDE("X")
2231              else:              else:
2232                  Xout=X[:]                  Xout=X[:]
2233              Xi=self.__getXi()              Xi=self.__getXi()
# Line 2228  class AdvectivePDE(LinearPDE): Line 2246  class AdvectivePDE(LinearPDE):
2246                               Xout[i,j]+=tmp*C[p,i,j]                               Xout[i,j]+=tmp*C[p,i,j]
2247                               # Xout=X_out+Xi*util.inner(Y,C,offset=1)                               # Xout=X_out+Xi*util.inner(Y,C,offset=1)
2248              else:              else:
2249                   tmp=Xi*Y                if not C.isEmpty() and not B.isEmpty():
2250                   for j in range(self.getDim()):                  Xout+=(Xi*Y)*(C-B)
2251                      if not C.isEmpty() and not B.isEmpty():                elif C.isEmpty():
2252                         Xout[j]+=tmp*(C[j]-B[j])                  Xout-=(Xi*Y)*B
2253                         # Xout=Xout+Xi*Y*(C-B)                else:
2254                      elif C.isEmpty():                  Xout+=(Xi*Y)*C
                        Xout[j]-=tmp*B[j]  
                        # Xout=Xout-Xi*Y*B  
                     else:  
                        Xout[j]+=tmp*C[j]  
                        # Xout=Xout+Xi*Y*C  
2255           return Xout           return Xout
2256       elif name == "Y" :       elif name == "Y" :
2257           return self.getCoefficient("Y")           return self.getCoefficient("Y")
# Line 2257  class AdvectivePDE(LinearPDE): Line 2270  class AdvectivePDE(LinearPDE):
2270       else:       else:
2271          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
2272    
 class AdvectionDiffusion(LinearPDE):  
    """  
    Class to define PDE equation of the unisotropic advection-diffusion problem, which is genear L{LinearPDE} of the form  
   
    M{S{omega}*u + inner(v,grad(u))- grad(matrixmult(k_bar,grad(u))[j])[j] = f}  
   
    with natural boundary conditons  
   
    M{n[j]*matrixmult(k,grad(u))[j] = g- S{alpha}u }  
   
    and constraints:  
   
    M{u=r} where M{q>0}  
   
    and  
   
    M{k_bar[i,j]=k[i,j]+upwind[i]*upwind[j]}  
   
    """  
   
    def __init__(self,domain,debug=False):  
      """  
      initializes a new Poisson equation  
   
      @param domain: domain of the PDE  
      @type domain: L{Domain<escript.Domain>}  
      @param debug: if True debug informations are printed.  
   
      """  
      super(AdvectionDiffusion, self).__init__(domain,1,1,debug)  
      self.COEFFICIENTS={"omega": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),  
                         "k": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_DIM,PDECoefficient.BY_DIM),PDECoefficient.OPERATOR),  
                         "f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),  
                         "v": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_DIM,),PDECoefficient.OPERATOR),  
                         "upwind": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_DIM,),PDECoefficient.OPERATOR),  
                         "alpha": PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),  
                         "g": PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),  
                         "r": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH),  
                         "q": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}  
   
    def setValue(self,**coefficients):  
      """  
      sets new values to coefficients  
   
      @param coefficients: new values assigned to coefficients  
      @keyword omega: value for coefficient M{S{omega}}  
      @type omega: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.  
      @keyword k: value for coefficient M{k}  
      @type k: any type that can be casted to L{Tensor<escript.Tensor>} object on L{Function<escript.Function>}.  
      @keyword v: value for coefficient M{v}  
      @type v: any type that can be casted to L{Vector<escript.Vector>} object on L{Function<escript.Function>}.  
      @keyword upwind: value for upwind term M{upwind}  
      @type upwind: any type that can be casted to L{Vector<escript.Vector>} object on L{Function<escript.Function>}.  
      @keyword f: value for right hand side M{f}  
      @type f: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.  
      @keyword alpha: value for right hand side M{S{alpha}}  
      @type alpha: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.  
      @keyword g: value for right hand side M{g}  
      @type g: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.  
      @keyword r: prescribed values M{r} for the solution in constraints.  
      @type r: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}  
                depending of reduced order is used for the representation of the equation.  
      @keyword q: mask for location of constraints  
      @type q: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}  
                depending of reduced order is used for the representation of the equation.  
      @raise IllegalCoefficient: if an unknown coefficient keyword is used.  
      """  
      super(AdvectionDiffusion, self).setValue(**coefficients)  
   
    def getCoefficientOfGeneralPDE(self,name):  
      """  
      return the value of the coefficient name of the general PDE  
   
      @param name: name of the coefficient requested.  
      @type name: C{string}  
      @return: the value of the coefficient  name  
      @rtype: L{Data<escript.Data>}  
      @raise IllegalCoefficient: if name is not one of coefficients  
                   "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}.  
      @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.  
      """  
      if name == "A" :  
          return self.getCoefficient("k")+util.outer(self.getCoefficient("upwind"),self.getCoefficient("upwind"))  
      elif name == "B" :  
          return escript.Data()  
      elif name == "C" :  
          return self.getCoefficient("v")  
      elif name == "D" :  
          return self.getCoefficient("omega")  
      elif name == "X" :  
          return escript.Data()  
      elif name == "Y" :  
          return self.getCoefficient("f")  
      elif name == "d" :  
          return self.getCoefficient("alpha")  
      elif name == "y" :  
          return self.getCoefficient("g")  
      elif name == "d_contact" :  
          return escript.Data()  
      elif name == "y_contact" :  
          return escript.Data()  
      elif name == "r" :  
          return self.getCoefficient("r")  
      elif name == "q" :  
          return self.getCoefficient("q")  
      else:  
         raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name  
   
   
2273  # $Log$  # $Log$
2274  # Revision 1.14  2005/09/22 01:54:57  jgs  # Revision 1.14  2005/09/22 01:54:57  jgs
2275  # Merge of development branch dev-02 back to main trunk on 2005-09-22  # Merge of development branch dev-02 back to main trunk on 2005-09-22

Legend:
Removed from v.410  
changed lines
  Added in v.969

  ViewVC Help
Powered by ViewVC 1.1.26