/[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 425 by gross, Tue Jan 10 04:10:39 2006 UTC revision 791 by bcumming, Thu Jul 27 00:37:57 2006 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 419  class LinearPDE(object): Line 414  class LinearPDE(object):
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 ITERATIVE: The default iterative solver     @cvar ITERATIVE: The default iterative solver
417       @cvar AMG: algebraic multi grid
418       @cvar RILU: recursive ILU
419    
420     """     """
421     DEFAULT= 0     DEFAULT= 0
# Line 443  class LinearPDE(object): Line 440  class LinearPDE(object):
440     UMFPACK= 16     UMFPACK= 16
441     ITERATIVE= 20     ITERATIVE= 20
442     PASO= 21     PASO= 21
443       AMG= 22
444       RILU = 23
445    
446     __TOL=1.e-13     SMALL_TOLERANCE=1.e-13
447     __PACKAGE_KEY="package"     __PACKAGE_KEY="package"
448     __METHOD_KEY="method"     __METHOD_KEY="method"
449     __SYMMETRY_KEY="symmetric"     __SYMMETRY_KEY="symmetric"
# Line 665  class LinearPDE(object): Line 664  class LinearPDE(object):
664       @rtype: L{Data<escript.Data>}       @rtype: L{Data<escript.Data>}
665       """       """
666       if u==None:       if u==None:
667            return self.getOperator()*self.getSolution()          return self.getOperator()*self.getSolution()
668       else:       else:
669          self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())          return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())
670    
671     def getResidual(self,u=None):     def getResidual(self,u=None):
672       """       """
# Line 699  class LinearPDE(object): Line 698  class LinearPDE(object):
698        else:        else:
699           A=self.getCoefficientOfGeneralPDE("A")           A=self.getCoefficientOfGeneralPDE("A")
700           if not A.isEmpty():           if not A.isEmpty():
701              tol=util.Lsup(A)*self.__TOL              tol=util.Lsup(A)*self.SMALL_TOLERANCE
702              if self.getNumSolutions()>1:              if self.getNumSolutions()>1:
703                 for i in range(self.getNumEquations()):                 for i in range(self.getNumEquations()):
704                    for j in range(self.getDim()):                    for j in range(self.getDim()):
# Line 723  class LinearPDE(object): Line 722  class LinearPDE(object):
722              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"
723              out=False              out=False
724           elif not B.isEmpty() and not C.isEmpty():           elif not B.isEmpty() and not C.isEmpty():
725              tol=(util.Lsup(B)+util.Lsup(C))*self.__TOL/2.              tol=(util.Lsup(B)+util.Lsup(C))*self.SMALL_TOLERANCE/2.
726              if self.getNumSolutions()>1:              if self.getNumSolutions()>1:
727                 for i in range(self.getNumEquations()):                 for i in range(self.getNumEquations()):
728                     for j in range(self.getDim()):                     for j in range(self.getDim()):
# Line 739  class LinearPDE(object): Line 738  class LinearPDE(object):
738           if self.getNumSolutions()>1:           if self.getNumSolutions()>1:
739             D=self.getCoefficientOfGeneralPDE("D")             D=self.getCoefficientOfGeneralPDE("D")
740             if not D.isEmpty():             if not D.isEmpty():
741               tol=util.Lsup(D)*self.__TOL               tol=util.Lsup(D)*self.SMALL_TOLERANCE
742               for i in range(self.getNumEquations()):               for i in range(self.getNumEquations()):
743                  for k in range(self.getNumSolutions()):                  for k in range(self.getNumSolutions()):
744                    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 746  class LinearPDE(object):
746                        out=False                        out=False
747             d=self.getCoefficientOfGeneralPDE("d")             d=self.getCoefficientOfGeneralPDE("d")
748             if not d.isEmpty():             if not d.isEmpty():
749               tol=util.Lsup(d)*self.__TOL               tol=util.Lsup(d)*self.SMALL_TOLERANCE
750               for i in range(self.getNumEquations()):               for i in range(self.getNumEquations()):
751                  for k in range(self.getNumSolutions()):                  for k in range(self.getNumSolutions()):
752                    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 754  class LinearPDE(object):
754                        out=False                        out=False
755             d_contact=self.getCoefficientOfGeneralPDE("d_contact")             d_contact=self.getCoefficientOfGeneralPDE("d_contact")
756             if not d_contact.isEmpty():             if not d_contact.isEmpty():
757               tol=util.Lsup(d_contact)*self.__TOL               tol=util.Lsup(d_contact)*self.SMALL_TOLERANCE
758               for i in range(self.getNumEquations()):               for i in range(self.getNumEquations()):
759                  for k in range(self.getNumSolutions()):                  for k in range(self.getNumSolutions()):
760                    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 820  class LinearPDE(object):
820         sets a new solver         sets a new solver
821    
822         @param solver: sets a new solver method.         @param solver: sets a new solver method.
823         @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}
824         @param preconditioner: sets a new solver method.         @param preconditioner: sets a new solver method.
825         @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}
826         """         """
827         if solver==None: solve=self.DEFAULT         if solver==None: solve=self.DEFAULT
828         if preconditioner==None: preconditioner=self.DEFAULT         if preconditioner==None: preconditioner=self.DEFAULT
# Line 856  class LinearPDE(object): Line 855  class LinearPDE(object):
855         elif m[0]==self.GMRES: method= "GMRES"         elif m[0]==self.GMRES: method= "GMRES"
856         elif m[0]==self.PRES20: method= "PRES20"         elif m[0]==self.PRES20: method= "PRES20"
857         elif m[0]==self.LUMPING: method= "LUMPING"         elif m[0]==self.LUMPING: method= "LUMPING"
858           elif m[0]==self.AMG: method= "AMG"
859         if m[1]==self.DEFAULT: method+="+DEFAULT"         if m[1]==self.DEFAULT: method+="+DEFAULT"
860         elif m[1]==self.JACOBI: method+= "+JACOBI"         elif m[1]==self.JACOBI: method+= "+JACOBI"
861         elif m[1]==self.ILU0: method+= "+ILU0"         elif m[1]==self.ILU0: method+= "+ILU0"
862         elif m[1]==self.ILUT: method+= "+ILUT"         elif m[1]==self.ILUT: method+= "+ILUT"
863         elif m[1]==self.SSOR: method+= "+SSOR"         elif m[1]==self.SSOR: method+= "+SSOR"
864           elif m[1]==self.AMG: method+= "+AMG"
865           elif m[1]==self.RILU: method+= "+RILU"
866         if p==self.DEFAULT: package="DEFAULT"         if p==self.DEFAULT: package="DEFAULT"
867         elif p==self.PASO: package= "PASO"         elif p==self.PASO: package= "PASO"
868         elif p==self.MKL: package= "MKL"         elif p==self.MKL: package= "MKL"
# Line 883  class LinearPDE(object): Line 885  class LinearPDE(object):
885         """         """
886         sets a new solver package         sets a new solver package
887    
888         @param solver: sets a new solver method.         @param package: sets a new solver method.
889         @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}
890         """         """
891         if package==None: package=self.DEFAULT         if package==None: package=self.DEFAULT
892         if not package==self.getSolverPackage():         if not package==self.getSolverPackage():
893             self.__solver_method=solver             self.__solver_package=package
894             self.__checkMatrixType()             self.__checkMatrixType()
895             self.trace("New solver is %s"%self.getSolverMethodName())             self.trace("New solver is %s"%self.getSolverMethodName())
896    
# Line 1511  class LinearPDE(object): Line 1513  class LinearPDE(object):
1513           r=self.getCoefficientOfGeneralPDE("r")           r=self.getCoefficientOfGeneralPDE("r")
1514           homogeneous_constraint=True           homogeneous_constraint=True
1515           if not q.isEmpty() and not r.isEmpty():           if not q.isEmpty() and not r.isEmpty():
1516               if util.Lsup(q*r)>=1.e-13*util.Lsup(r):               if util.Lsup(q*r)>0.:
1517                 self.trace("Inhomogeneous constraint detected.")                 self.trace("Inhomogeneous constraint detected.")
1518                 self.__invalidateSystem()                 self.__invalidateSystem()
1519    
# Line 1525  class LinearPDE(object): Line 1527  class LinearPDE(object):
1527         if not self.__operator_is_Valid or not self.__righthandside_isValid:         if not self.__operator_is_Valid or not self.__righthandside_isValid:
1528            if self.isUsingLumping():            if self.isUsingLumping():
1529                if not self.__operator_is_Valid:                if not self.__operator_is_Valid:
1530                   if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution(): raise TypeError,"Lumped matrix requires same order for equations and unknowns"                   if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():
1531                   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"
1532                   if not self.getCoefficientOfGeneralPDE("B").isEmpty(): raise Warning,"Using coefficient B in lumped matrix can produce wrong results"                   if not self.getCoefficientOfGeneralPDE("A").isEmpty():
1533                   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."
1534                   mat=self.__getNewOperator()                   if not self.getCoefficientOfGeneralPDE("B").isEmpty():
1535                   self.getDomain().addPDEToSystem(mat,escript.Data(), \                        raise ValueError,"coefficient A in lumped matrix may not be present."
1536                             self.getCoefficientOfGeneralPDE("A"), \                   if not self.getCoefficientOfGeneralPDE("C").isEmpty():
1537                             self.getCoefficientOfGeneralPDE("B"), \                        raise ValueError,"coefficient A in lumped matrix may not be present."
1538                             self.getCoefficientOfGeneralPDE("C"), \                   D=self.getCoefficientOfGeneralPDE("D")
1539                             self.getCoefficientOfGeneralPDE("D"), \                   if not D.isEmpty():
1540                             escript.Data(), \                       if self.getNumSolutions()>1:
1541                             escript.Data(), \                          #D_times_e=util.matrixmult(D,numarray.ones((self.getNumSolutions(),)))
1542                             self.getCoefficientOfGeneralPDE("d"), \                          D_times_e=util.matrix_mult(D,numarray.ones((self.getNumSolutions(),)))
1543                             escript.Data(),\                       else:
1544                             self.getCoefficientOfGeneralPDE("d_contact"), \                          D_times_e=D
1545                             escript.Data())                   else:
1546                   self.__operator=1./(mat*escript.Data(1,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True))                      D_times_e=escript.Data()
1547                   del mat                   d=self.getCoefficientOfGeneralPDE("d")
1548                     if not d.isEmpty():
1549                         if self.getNumSolutions()>1:
1550                            #d_times_e=util.matrixmult(d,numarray.ones((self.getNumSolutions(),)))
1551                            d_times_e=util.matrix_mult(d,numarray.ones((self.getNumSolutions(),)))
1552                         else:
1553                            d_times_e=d
1554                     else:
1555                        d_times_e=escript.Data()
1556                     d_contact=self.getCoefficientOfGeneralPDE("d_contact")
1557                     if not d_contact.isEmpty():
1558                         if self.getNumSolutions()>1:
1559                            d_contact_times_e=util.matrixmult(d_contact,numarray.ones((self.getNumSolutions(),)))
1560                         else:
1561                            d_contact_times_e=d_contact
1562                     else:
1563                        d_contact_times_e=escript.Data()
1564        
1565                     self.__operator=self.__getNewRightHandSide()
1566                     self.getDomain().addPDEToRHS(self.__operator, \
1567                                                  escript.Data(), \
1568                                                  D_times_e, \
1569                                                  d_times_e,\
1570                                                  d_contact_times_e)
1571                     self.__operator=1./self.__operator
1572                   self.trace("New lumped operator has been built.")                   self.trace("New lumped operator has been built.")
1573                   self.__operator_is_Valid=True                   self.__operator_is_Valid=True
1574                if not self.__righthandside_isValid:                if not self.__righthandside_isValid:
# Line 1783  class LameEquation(LinearPDE): Line 1809  class LameEquation(LinearPDE):
1809     """     """
1810     Class to define a Lame equation problem:     Class to define a Lame equation problem:
1811    
1812     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] }
1813    
1814     with natural boundary conditons:     with natural boundary conditons:
1815    
1816     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] }
1817    
1818     and constraints:     and constraints:
1819    
# Line 1807  class LameEquation(LinearPDE): Line 1833  class LameEquation(LinearPDE):
1833                            "q"            : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}                            "q"            : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}
1834        self.setSymmetryOn()        self.setSymmetryOn()
1835    
1836     def setValue(self,**coefficients):     def setValues(self,**coefficients):
1837       """       """
1838       sets new values to coefficients       sets new values to coefficients
1839    
# Line 1830  class LameEquation(LinearPDE): Line 1856  class LameEquation(LinearPDE):
1856                 depending of reduced order is used for the representation of the equation.                 depending of reduced order is used for the representation of the equation.
1857       @raise IllegalCoefficient: if an unknown coefficient keyword is used.       @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1858       """       """
1859       super(LameEquation, self).setValue(**coefficients)       super(LameEquation, self).setValues(**coefficients)
1860    
1861     def getCoefficientOfGeneralPDE(self,name):     def getCoefficientOfGeneralPDE(self,name):
1862       """       """
# Line 1944  class AdvectivePDE(LinearPDE): Line 1970  class AdvectivePDE(LinearPDE):
1970           self.__xi=xi           self.__xi=xi
1971        self.__Xi=escript.Data()        self.__Xi=escript.Data()
1972    
1973     def setValue(**coefficients):     def setValue(self,**coefficients):
1974        """        """
1975        sets new values to coefficients        sets new values to coefficients
1976    
# Line 2028  class AdvectivePDE(LinearPDE): Line 2054  class AdvectivePDE(LinearPDE):
2054       """       """
2055       return escript.Scalar(0.5,P.getFunctionSpace())       return escript.Scalar(0.5,P.getFunctionSpace())
2056    
    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.  
   
2057     def __getXi(self):     def __getXi(self):
2058        if self.__Xi.isEmpty():        if self.__Xi.isEmpty():
2059           B=self.getCoefficient("B")           B=self.getCoefficient("B")
# Line 2044  class AdvectivePDE(LinearPDE): Line 2063  class AdvectivePDE(LinearPDE):
2063           self.__Xi=escript.Scalar(0.,self.getFunctionSpaceForCoefficient("A"))           self.__Xi=escript.Scalar(0.,self.getFunctionSpaceForCoefficient("A"))
2064           if not C.isEmpty() or not B.isEmpty():           if not C.isEmpty() or not B.isEmpty():
2065              if not C.isEmpty() and not B.isEmpty():              if not C.isEmpty() and not B.isEmpty():
                 flux2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))  
2066                  if self.getNumEquations()>1:                  if self.getNumEquations()>1:
2067                     if self.getNumSolutions()>1:                     if self.getNumSolutions()>1:
2068                          flux2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))
2069                        for i in range(self.getNumEquations()):                        for i in range(self.getNumEquations()):
2070                           for k in range(self.getNumSolutions()):                           for k in range(self.getNumSolutions()):
2071                              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
2072                          length_of_flux=util.sqrt(flux2)
2073                        # flux=C-util.reorderComponents(B,[0,2,1])                        # flux=C-util.reorderComponents(B,[0,2,1])
2074                     else:                     else:
2075                          flux2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))
2076                        for i in range(self.getNumEquations()):                        for i in range(self.getNumEquations()):
2077                           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
2078                          length_of_flux=util.sqrt(flux2)
2079                        # flux=C-B                        # flux=C-B
2080                  else:                  else:
2081                     if self.getNumSolutions()>1:                     if self.getNumSolutions()>1:
2082                          flux2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))
2083                        for k in range(self.getNumSolutions()):                        for k in range(self.getNumSolutions()):
2084                           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
2085                        # flux=C-util.reorderComponents(B,[1,0])                        # flux=C-util.reorderComponents(B,[1,0])
2086                          length_of_flux=util.sqrt(flux2)
2087                     else:                     else:
2088                        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)  
2089              elif C.isEmpty():              elif C.isEmpty():
2090                length_of_flux=util.length(B)                length_of_flux=util.length(B)
               #flux=B  
2091              else:              else:
2092                length_of_flux=util.length(C)                length_of_flux=util.length(C)
               #flux=C  
   
             #length_of_flux=util.length(flux)  
2093              flux_max=util.Lsup(length_of_flux)              flux_max=util.Lsup(length_of_flux)
2094              if flux_max>0.:              if flux_max>0.:
2095                 # length_of_A=util.inner(flux,util.tensormutiply(A,flux))                if A.isEmpty():
2096                 length_of_A=util.length(A)                    inv_A=1./self.SMALL_TOLERANCE
2097                 A_max=util.Lsup(length_of_A)                    peclet_number=escript.Scalar(inv_A,length_of_flux.getFunctionSpace())
2098                 if A_max>0:                    xi=self.__xi(self,peclet_number)
2099                      inv_A=1./(length_of_A+A_max*self.__TOL)                else:
2100                 else:                    # length_of_A=util.inner(flux,util.tensormutiply(A,flux))
2101                      inv_A=1./self.__TOL                    length_of_A=util.length(A)
2102                 peclet_number=length_of_flux*h/2*inv_A                    A_max=util.Lsup(length_of_A)
2103                 xi=self.__xi(peclet_number)                    if A_max>0:
2104                 self.__Xi=h*xi/(length_of_flux+flux_max*self.__TOL)                         inv_A=1./(length_of_A+A_max*self.SMALL_TOLERANCE)
2105                 self.trace("preclet number = %e"%util.Lsup(peclet_number))                    else:
2106                           inv_A=1./self.SMALL_TOLERANCE
2107                      peclet_number=length_of_flux*h/2*inv_A
2108                      xi=self.__xi(self,peclet_number)
2109                  self.__Xi=h*xi/(length_of_flux+flux_max*self.SMALL_TOLERANCE)
2110                  self.trace("preclet number = %e"%util.Lsup(peclet_number))
2111                else:
2112                  self.__Xi=escript.Scalar(0.,length_of_flux.getFunctionSpace())
2113        return self.__Xi        return self.__Xi
2114    
2115    
# Line 2111  class AdvectivePDE(LinearPDE): Line 2136  class AdvectivePDE(LinearPDE):
2136              Aout=A              Aout=A
2137           else:           else:
2138              if A.isEmpty():              if A.isEmpty():
2139                 Aout=self.createNewCoefficient("A")                 Aout=self.createCoefficientOfGeneralPDE("A")
2140              else:              else:
2141                 Aout=A[:]                 Aout=A[:]
2142              Xi=self.__getXi()              Xi=self.__getXi()
# Line 2131  class AdvectivePDE(LinearPDE): Line 2156  class AdvectivePDE(LinearPDE):
2156                                 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]
2157                                 # 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)
2158              else:              else:
2159                  for j in range(self.getDim()):                 if not C.isEmpty() and not B.isEmpty():
2160                     for l in range(self.getDim()):                     delta=(C-B)
2161                        if not C.isEmpty() and not B.isEmpty():                     Aout+=util.outer(Xi*delta,delta)
2162                            Aout[j,l]+=Xi*(C[j]-B[j])*(C[l]-B[l])                 elif not B.isEmpty():
2163                        elif C.isEmpty():                     Aout+=util.outer(Xi*B,B)
2164                            Aout[j,l]+=Xi*B[j]*B[l]                 elif not C.isEmpty():
2165                        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)  
2166           return Aout           return Aout
2167       elif name == "B" :       elif name == "B" :
2168             # return self.getCoefficient("B")
2169           B=self.getCoefficient("B")           B=self.getCoefficient("B")
2170           C=self.getCoefficient("C")           C=self.getCoefficient("C")
2171           D=self.getCoefficient("D")           D=self.getCoefficient("D")
# Line 2156  class AdvectivePDE(LinearPDE): Line 2174  class AdvectivePDE(LinearPDE):
2174           else:           else:
2175              Xi=self.__getXi()              Xi=self.__getXi()
2176              if B.isEmpty():              if B.isEmpty():
2177                  Bout=self.createNewCoefficient("B")                  Bout=self.createCoefficientOfGeneralPDE("B")
2178              else:              else:
2179                  Bout=B[:]                  Bout=B[:]
2180              if self.getNumEquations()>1:              if self.getNumEquations()>1:
# Line 2168  class AdvectivePDE(LinearPDE): Line 2186  class AdvectivePDE(LinearPDE):
2186                             Bout[i,j,k]+=tmp*C[p,i,j]                             Bout[i,j,k]+=tmp*C[p,i,j]
2187                             # 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)
2188              else:              else:
2189                 tmp=Xi*D                 Bout+=(Xi*D)*C
                for j in range(self.getDim()): Bout[j]+=tmp*C[j]  
                # Bout=Bout+Xi*D*C  
2190           return Bout           return Bout
2191       elif name == "C" :       elif name == "C" :
2192             # return self.getCoefficient("C")
2193           B=self.getCoefficient("B")           B=self.getCoefficient("B")
2194           C=self.getCoefficient("C")           C=self.getCoefficient("C")
2195           D=self.getCoefficient("D")           D=self.getCoefficient("D")
# Line 2181  class AdvectivePDE(LinearPDE): Line 2198  class AdvectivePDE(LinearPDE):
2198           else:           else:
2199              Xi=self.__getXi()              Xi=self.__getXi()
2200              if C.isEmpty():              if C.isEmpty():
2201                  Cout=self.createNewCoefficient("C")                  Cout=self.createCoefficientOfGeneralPDE("C")
2202              else:              else:
2203                  Cout=C[:]                  Cout=C[:]
2204              if self.getNumEquations()>1:              if self.getNumEquations()>1:
# Line 2193  class AdvectivePDE(LinearPDE): Line 2210  class AdvectivePDE(LinearPDE):
2210                                   Cout[i,k,l]+=tmp*B[p,l,i]                                   Cout[i,k,l]+=tmp*B[p,l,i]
2211                                   # Cout=Cout+Xi*B[p,l,i]*D[p,k]                                   # Cout=Cout+Xi*B[p,l,i]*D[p,k]
2212              else:              else:
2213                 tmp=Xi*D                 Cout+=(Xi*D)*B
                for j in range(self.getDim()): Cout[j]+=tmp*B[j]  
                # Cout=Cout+tmp*D*B  
2214           return Cout           return Cout
2215       elif name == "D" :       elif name == "D" :
2216           return self.getCoefficient("D")           return self.getCoefficient("D")
2217       elif name == "X" :       elif name == "X" :
2218             # return self.getCoefficient("X")
2219           X=self.getCoefficient("X")           X=self.getCoefficient("X")
2220           Y=self.getCoefficient("Y")           Y=self.getCoefficient("Y")
2221           B=self.getCoefficient("B")           B=self.getCoefficient("B")
# Line 2208  class AdvectivePDE(LinearPDE): Line 2224  class AdvectivePDE(LinearPDE):
2224              Xout=X              Xout=X
2225           else:           else:
2226              if X.isEmpty():              if X.isEmpty():
2227                  Xout=self.createNewCoefficient("X")                  Xout=self.createCoefficientOfGeneralPDE("X")
2228              else:              else:
2229                  Xout=X[:]                  Xout=X[:]
2230              Xi=self.__getXi()              Xi=self.__getXi()
# Line 2227  class AdvectivePDE(LinearPDE): Line 2243  class AdvectivePDE(LinearPDE):
2243                               Xout[i,j]+=tmp*C[p,i,j]                               Xout[i,j]+=tmp*C[p,i,j]
2244                               # Xout=X_out+Xi*util.inner(Y,C,offset=1)                               # Xout=X_out+Xi*util.inner(Y,C,offset=1)
2245              else:              else:
2246                   tmp=Xi*Y                if not C.isEmpty() and not B.isEmpty():
2247                   for j in range(self.getDim()):                  Xout+=(Xi*Y)*(C-B)
2248                      if not C.isEmpty() and not B.isEmpty():                elif C.isEmpty():
2249                         Xout[j]+=tmp*(C[j]-B[j])                  Xout-=(Xi*Y)*B
2250                         # Xout=Xout+Xi*Y*(C-B)                else:
2251                      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  
2252           return Xout           return Xout
2253       elif name == "Y" :       elif name == "Y" :
2254           return self.getCoefficient("Y")           return self.getCoefficient("Y")
# Line 2256  class AdvectivePDE(LinearPDE): Line 2267  class AdvectivePDE(LinearPDE):
2267       else:       else:
2268          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
2269    
 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  
   
   
2270  # $Log$  # $Log$
2271  # Revision 1.14  2005/09/22 01:54:57  jgs  # Revision 1.14  2005/09/22 01:54:57  jgs
2272  # 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.425  
changed lines
  Added in v.791

  ViewVC Help
Powered by ViewVC 1.1.26