/[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 345 by gross, Tue Dec 13 05:23:45 2005 UTC revision 614 by elspeth, Wed Mar 22 01:37:07 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 22  classs which are also derived form the L Line 12  classs which are also derived form the L
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 __license__: licence agreement
16  @var __url__: url entry point on documentation  @var __url__: url entry point on documentation
17  @var __version__: version  @var __version__: version
18  @var __date__: date of the version  @var __date__: date of the version
# Line 33  import util Line 23  import util
23  import numarray  import numarray
24    
25  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
26  __licence__="contact: esys@access.uq.edu.au"  __copyright__="""  Copyright (c) 2006 by ACcESS MNRF
27                        http://www.access.edu.au
28                    Primary Business: Queensland, Australia"""
29    __license__="""Licensed under the Open Software License version 3.0
30                 http://www.opensource.org/licenses/osl-3.0.php"""
31  __url__="http://www.iservo.edu.au/esys/escript"  __url__="http://www.iservo.edu.au/esys/escript"
32  __version__="$Revision$"  __version__="$Revision$"
33  __date__="$Date$"  __date__="$Date$"
# Line 419  class LinearPDE(object): Line 413  class LinearPDE(object):
413     @cvar MKL: Intel's MKL solver library     @cvar MKL: Intel's MKL solver library
414     @cvar UMFPACK: the UMFPACK library     @cvar UMFPACK: the UMFPACK library
415     @cvar ITERATIVE: The default iterative solver     @cvar ITERATIVE: The default iterative solver
416       @cvar AMG: algebraic multi grid
417       @cvar RILU: recursive ILU
418    
419     """     """
420     DEFAULT= 0     DEFAULT= 0
# Line 443  class LinearPDE(object): Line 439  class LinearPDE(object):
439     UMFPACK= 16     UMFPACK= 16
440     ITERATIVE= 20     ITERATIVE= 20
441     PASO= 21     PASO= 21
442       AMG= 22
443       RILU = 23
444    
445     __TOL=1.e-13     SMALL_TOLERANCE=1.e-13
446     __PACKAGE_KEY="package"     __PACKAGE_KEY="package"
447     __METHOD_KEY="method"     __METHOD_KEY="method"
448     __SYMMETRY_KEY="symmetric"     __SYMMETRY_KEY="symmetric"
449     __TOLERANCE_KEY="tolerance"     __TOLERANCE_KEY="tolerance"
450       __PRECONDITIONER_KEY="preconditioner"
451    
452    
453     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 497  class LinearPDE(object):
497       self.__tolerance=1.e-8       self.__tolerance=1.e-8
498       self.__solver_method=self.DEFAULT       self.__solver_method=self.DEFAULT
499       self.__solver_package=self.DEFAULT       self.__solver_package=self.DEFAULT
500         self.__preconditioner=self.DEFAULT
501       self.__matrix_type=self.__domain.getSystemMatrixTypeId(self.DEFAULT,self.DEFAULT,False)       self.__matrix_type=self.__domain.getSystemMatrixTypeId(self.DEFAULT,self.DEFAULT,False)
502       self.__sym=False       self.__sym=False
503    
# Line 697  class LinearPDE(object): Line 697  class LinearPDE(object):
697        else:        else:
698           A=self.getCoefficientOfGeneralPDE("A")           A=self.getCoefficientOfGeneralPDE("A")
699           if not A.isEmpty():           if not A.isEmpty():
700              tol=util.Lsup(A)*self.__TOL              tol=util.Lsup(A)*self.SMALL_TOLERANCE
701              if self.getNumSolutions()>1:              if self.getNumSolutions()>1:
702                 for i in range(self.getNumEquations()):                 for i in range(self.getNumEquations()):
703                    for j in range(self.getDim()):                    for j in range(self.getDim()):
# Line 721  class LinearPDE(object): Line 721  class LinearPDE(object):
721              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"
722              out=False              out=False
723           elif not B.isEmpty() and not C.isEmpty():           elif not B.isEmpty() and not C.isEmpty():
724              tol=(util.Lsup(B)+util.Lsup(C))*self.__TOL/2.              tol=(util.Lsup(B)+util.Lsup(C))*self.SMALL_TOLERANCE/2.
725              if self.getNumSolutions()>1:              if self.getNumSolutions()>1:
726                 for i in range(self.getNumEquations()):                 for i in range(self.getNumEquations()):
727                     for j in range(self.getDim()):                     for j in range(self.getDim()):
# Line 737  class LinearPDE(object): Line 737  class LinearPDE(object):
737           if self.getNumSolutions()>1:           if self.getNumSolutions()>1:
738             D=self.getCoefficientOfGeneralPDE("D")             D=self.getCoefficientOfGeneralPDE("D")
739             if not D.isEmpty():             if not D.isEmpty():
740               tol=util.Lsup(D)*self.__TOL               tol=util.Lsup(D)*self.SMALL_TOLERANCE
741               for i in range(self.getNumEquations()):               for i in range(self.getNumEquations()):
742                  for k in range(self.getNumSolutions()):                  for k in range(self.getNumSolutions()):
743                    if util.Lsup(D[i,k]-D[k,i])>tol:                    if util.Lsup(D[i,k]-D[k,i])>tol:
# Line 745  class LinearPDE(object): Line 745  class LinearPDE(object):
745                        out=False                        out=False
746             d=self.getCoefficientOfGeneralPDE("d")             d=self.getCoefficientOfGeneralPDE("d")
747             if not d.isEmpty():             if not d.isEmpty():
748               tol=util.Lsup(d)*self.__TOL               tol=util.Lsup(d)*self.SMALL_TOLERANCE
749               for i in range(self.getNumEquations()):               for i in range(self.getNumEquations()):
750                  for k in range(self.getNumSolutions()):                  for k in range(self.getNumSolutions()):
751                    if util.Lsup(d[i,k]-d[k,i])>tol:                    if util.Lsup(d[i,k]-d[k,i])>tol:
# Line 753  class LinearPDE(object): Line 753  class LinearPDE(object):
753                        out=False                        out=False
754             d_contact=self.getCoefficientOfGeneralPDE("d_contact")             d_contact=self.getCoefficientOfGeneralPDE("d_contact")
755             if not d_contact.isEmpty():             if not d_contact.isEmpty():
756               tol=util.Lsup(d_contact)*self.__TOL               tol=util.Lsup(d_contact)*self.SMALL_TOLERANCE
757               for i in range(self.getNumEquations()):               for i in range(self.getNumEquations()):
758                  for k in range(self.getNumSolutions()):                  for k in range(self.getNumSolutions()):
759                    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 772  class LinearPDE(object): Line 772  class LinearPDE(object):
772         @type verbose: C{bool}         @type verbose: C{bool}
773         @keyword reordering: reordering scheme to be used during elimination. Allowed values are         @keyword reordering: reordering scheme to be used during elimination. Allowed values are
774                              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}  
775         @keyword iter_max: maximum number of iteration steps allowed.         @keyword iter_max: maximum number of iteration steps allowed.
776         @keyword drop_tolerance: threshold for drupping in L{ILUT}         @keyword drop_tolerance: threshold for drupping in L{ILUT}
777         @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 784  class LinearPDE(object):
784               self.__solution=self.copyConstraint(f*mat)               self.__solution=self.copyConstraint(f*mat)
785            else:            else:
786               options[self.__TOLERANCE_KEY]=self.getTolerance()               options[self.__TOLERANCE_KEY]=self.getTolerance()
787               options[self.__METHOD_KEY]=self.getSolverMethod()               options[self.__METHOD_KEY]=self.getSolverMethod()[0]
788                 options[self.__PRECONDITIONER_KEY]=self.getSolverMethod()[1]
789               options[self.__PACKAGE_KEY]=self.getSolverPackage()               options[self.__PACKAGE_KEY]=self.getSolverPackage()
790               options[self.__SYMMETRY_KEY]=self.isSymmetric()               options[self.__SYMMETRY_KEY]=self.isSymmetric()
791               self.trace("PDE is resolved.")               self.trace("PDE is resolved.")
# Line 815  class LinearPDE(object): Line 814  class LinearPDE(object):
814     # =============================================================================     # =============================================================================
815     #   solver settings:     #   solver settings:
816     # =============================================================================     # =============================================================================
817     def setSolverMethod(self,solver=None):     def setSolverMethod(self,solver=None,preconditioner=None):
818         """         """
819         sets a new solver         sets a new solver
820    
821         @param solver: sets a new solver method.         @param solver: sets a new solver method.
822         @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}
823           @param preconditioner: sets a new solver method.
824           @type preconditioner: one of L{DEFAULT}, L{JACOBI} L{ILU0}, L{ILUT},L{SSOR}, L{RILU}
825         """         """
826         if solver==None: solve=self.DEFAULT         if solver==None: solve=self.DEFAULT
827         if not solver==self.getSolverMethod():         if preconditioner==None: preconditioner=self.DEFAULT
828           if not (solver,preconditioner)==self.getSolverMethod():
829             self.__solver_method=solver             self.__solver_method=solver
830               self.__preconditioner=preconditioner
831             self.__checkMatrixType()             self.__checkMatrixType()
832             self.trace("New solver is %s"%self.getSolverMethodName())             self.trace("New solver is %s"%self.getSolverMethodName())
833    
# Line 838  class LinearPDE(object): Line 841  class LinearPDE(object):
841    
842         m=self.getSolverMethod()         m=self.getSolverMethod()
843         p=self.getSolverPackage()         p=self.getSolverPackage()
844         if m==self.DEFAULT: method="DEFAULT"         method=""
845         elif m==self.DIRECT: method= "DIRECT"         if m[0]==self.DEFAULT: method="DEFAULT"
846         elif m==self.ITERATIVE: method= "ITERATIVE"         elif m[0]==self.DIRECT: method= "DIRECT"
847         elif m==self.CHOLEVSKY: method= "CHOLEVSKY"         elif m[0]==self.ITERATIVE: method= "ITERATIVE"
848         elif m==self.PCG: method= "PCG"         elif m[0]==self.CHOLEVSKY: method= "CHOLEVSKY"
849         elif m==self.CR: method= "CR"         elif m[0]==self.PCG: method= "PCG"
850         elif m==self.CGS: method= "CGS"         elif m[0]==self.CR: method= "CR"
851         elif m==self.BICGSTAB: method= "BICGSTAB"         elif m[0]==self.CGS: method= "CGS"
852         elif m==self.SSOR: method= "SSOR"         elif m[0]==self.BICGSTAB: method= "BICGSTAB"
853         elif m==self.GMRES: method= "GMRES"         elif m[0]==self.SSOR: method= "SSOR"
854         elif m==self.PRES20: method= "PRES20"         elif m[0]==self.GMRES: method= "GMRES"
855         elif m==self.LUMPING: method= "LUMPING"         elif m[0]==self.PRES20: method= "PRES20"
856         else : method="unknown"         elif m[0]==self.LUMPING: method= "LUMPING"
857           elif m[0]==self.AMG: method= "AMG"
858           if m[1]==self.DEFAULT: method+="+DEFAULT"
859           elif m[1]==self.JACOBI: method+= "+JACOBI"
860           elif m[1]==self.ILU0: method+= "+ILU0"
861           elif m[1]==self.ILUT: method+= "+ILUT"
862           elif m[1]==self.SSOR: method+= "+SSOR"
863           elif m[1]==self.AMG: method+= "+AMG"
864           elif m[1]==self.RILU: method+= "+RILU"
865         if p==self.DEFAULT: package="DEFAULT"         if p==self.DEFAULT: package="DEFAULT"
866         elif p==self.PASO: package= "PASO"         elif p==self.PASO: package= "PASO"
867         elif p==self.MKL: package= "MKL"         elif p==self.MKL: package= "MKL"
# Line 867  class LinearPDE(object): Line 878  class LinearPDE(object):
878         @return: the solver method currently be used.         @return: the solver method currently be used.
879         @rtype: C{int}         @rtype: C{int}
880         """         """
881         return self.__solver_method         return self.__solver_method,self.__preconditioner
882    
883     def setSolverPackage(self,package=None):     def setSolverPackage(self,package=None):
884         """         """
# Line 898  class LinearPDE(object): Line 909  class LinearPDE(object):
909        @return: True is lumping is currently used a solver method.        @return: True is lumping is currently used a solver method.
910        @rtype: C{bool}        @rtype: C{bool}
911        """        """
912        return self.getSolverMethod()==self.LUMPING        return self.getSolverMethod()[0]==self.LUMPING
913    
914     def setTolerance(self,tol=1.e-8):     def setTolerance(self,tol=1.e-8):
915         """         """
# Line 1091  class LinearPDE(object): Line 1102  class LinearPDE(object):
1102       """       """
1103       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.
1104       """       """
1105       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())
1106       if not new_matrix_type==self.__matrix_type:       if not new_matrix_type==self.__matrix_type:
1107           self.trace("Matrix type is now %d."%new_matrix_type)           self.trace("Matrix type is now %d."%new_matrix_type)
1108           self.__matrix_type=new_matrix_type           self.__matrix_type=new_matrix_type
# Line 1515  class LinearPDE(object): Line 1526  class LinearPDE(object):
1526         if not self.__operator_is_Valid or not self.__righthandside_isValid:         if not self.__operator_is_Valid or not self.__righthandside_isValid:
1527            if self.isUsingLumping():            if self.isUsingLumping():
1528                if not self.__operator_is_Valid:                if not self.__operator_is_Valid:
1529                   if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution(): raise TypeError,"Lumped matrix requires same order for equations and unknowns"                   if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():
1530                   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"
1531                   if not self.getCoefficientOfGeneralPDE("B").isEmpty(): raise Warning,"Using coefficient B in lumped matrix can produce wrong results"                   if not self.getCoefficientOfGeneralPDE("A").isEmpty():
1532                   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."
1533                   mat=self.__getNewOperator()                   if not self.getCoefficientOfGeneralPDE("B").isEmpty():
1534                   self.getDomain().addPDEToSystem(mat,escript.Data(), \                        raise ValueError,"coefficient A in lumped matrix may not be present."
1535                             self.getCoefficientOfGeneralPDE("A"), \                   if not self.getCoefficientOfGeneralPDE("C").isEmpty():
1536                             self.getCoefficientOfGeneralPDE("B"), \                        raise ValueError,"coefficient A in lumped matrix may not be present."
1537                             self.getCoefficientOfGeneralPDE("C"), \                   D=self.getCoefficientOfGeneralPDE("D")
1538                             self.getCoefficientOfGeneralPDE("D"), \                   if not D.isEmpty():
1539                             escript.Data(), \                       if self.getNumSolutions()>1:
1540                             escript.Data(), \                          D_times_e=util.matrixmult(D,numarray.ones((self.getNumSolutions(),)))
1541                             self.getCoefficientOfGeneralPDE("d"), \                       else:
1542                             escript.Data(),\                          D_times_e=D
1543                             self.getCoefficientOfGeneralPDE("d_contact"), \                   else:
1544                             escript.Data())                      D_times_e=escript.Data()
1545                   self.__operator=1./(mat*escript.Data(1,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True))                   d=self.getCoefficientOfGeneralPDE("d")
1546                   del mat                   if not d.isEmpty():
1547                         if self.getNumSolutions()>1:
1548                            d_times_e=util.matrixmult(d,numarray.ones((self.getNumSolutions(),)))
1549                         else:
1550                            d_times_e=d
1551                     else:
1552                        d_times_e=escript.Data()
1553                     d_contact=self.getCoefficientOfGeneralPDE("d_contact")
1554                     if not d_contact.isEmpty():
1555                         if self.getNumSolutions()>1:
1556                            d_contact_times_e=util.matrixmult(d_contact,numarray.ones((self.getNumSolutions(),)))
1557                         else:
1558                            d_contact_times_e=d_contact
1559                     else:
1560                        d_contact_times_e=escript.Data()
1561        
1562                     self.__operator=self.__getNewRightHandSide()
1563                     self.getDomain().addPDEToRHS(self.__operator, \
1564                                                  escript.Data(), \
1565                                                  D_times_e, \
1566                                                  d_times_e,\
1567                                                  d_contact_times_e)
1568                     self.__operator=1./self.__operator
1569                   self.trace("New lumped operator has been built.")                   self.trace("New lumped operator has been built.")
1570                   self.__operator_is_Valid=True                   self.__operator_is_Valid=True
1571                if not self.__righthandside_isValid:                if not self.__righthandside_isValid:
# Line 1797  class LameEquation(LinearPDE): Line 1830  class LameEquation(LinearPDE):
1830                            "q"            : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}                            "q"            : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}
1831        self.setSymmetryOn()        self.setSymmetryOn()
1832    
1833     def setValue(self,**coefficients):     def setValues(self,**coefficients):
1834       """       """
1835       sets new values to coefficients       sets new values to coefficients
1836    
# Line 1820  class LameEquation(LinearPDE): Line 1853  class LameEquation(LinearPDE):
1853                 depending of reduced order is used for the representation of the equation.                 depending of reduced order is used for the representation of the equation.
1854       @raise IllegalCoefficient: if an unknown coefficient keyword is used.       @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1855       """       """
1856       super(LameEquation, self).setValue(**coefficients)       super(LameEquation, self).setValues(**coefficients)
1857    
1858     def getCoefficientOfGeneralPDE(self,name):     def getCoefficientOfGeneralPDE(self,name):
1859       """       """
# Line 1934  class AdvectivePDE(LinearPDE): Line 1967  class AdvectivePDE(LinearPDE):
1967           self.__xi=xi           self.__xi=xi
1968        self.__Xi=escript.Data()        self.__Xi=escript.Data()
1969    
1970     def setValue(**coefficients):     def setValue(self,**coefficients):
1971        """        """
1972        sets new values to coefficients        sets new values to coefficients
1973    
# Line 2018  class AdvectivePDE(LinearPDE): Line 2051  class AdvectivePDE(LinearPDE):
2051       """       """
2052       return escript.Scalar(0.5,P.getFunctionSpace())       return escript.Scalar(0.5,P.getFunctionSpace())
2053    
    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.  
   
2054     def __getXi(self):     def __getXi(self):
2055        if self.__Xi.isEmpty():        if self.__Xi.isEmpty():
2056           B=self.getCoefficient("B")           B=self.getCoefficient("B")
# Line 2034  class AdvectivePDE(LinearPDE): Line 2060  class AdvectivePDE(LinearPDE):
2060           self.__Xi=escript.Scalar(0.,self.getFunctionSpaceForCoefficient("A"))           self.__Xi=escript.Scalar(0.,self.getFunctionSpaceForCoefficient("A"))
2061           if not C.isEmpty() or not B.isEmpty():           if not C.isEmpty() or not B.isEmpty():
2062              if not C.isEmpty() and not B.isEmpty():              if not C.isEmpty() and not B.isEmpty():
                 flux2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))  
2063                  if self.getNumEquations()>1:                  if self.getNumEquations()>1:
2064                     if self.getNumSolutions()>1:                     if self.getNumSolutions()>1:
2065                          flux2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))
2066                        for i in range(self.getNumEquations()):                        for i in range(self.getNumEquations()):
2067                           for k in range(self.getNumSolutions()):                           for k in range(self.getNumSolutions()):
2068                              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
2069                          length_of_flux=util.sqrt(flux2)
2070                        # flux=C-util.reorderComponents(B,[0,2,1])                        # flux=C-util.reorderComponents(B,[0,2,1])
2071                     else:                     else:
2072                          flux2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))
2073                        for i in range(self.getNumEquations()):                        for i in range(self.getNumEquations()):
2074                           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
2075                          length_of_flux=util.sqrt(flux2)
2076                        # flux=C-B                        # flux=C-B
2077                  else:                  else:
2078                     if self.getNumSolutions()>1:                     if self.getNumSolutions()>1:
2079                          flux2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))
2080                        for k in range(self.getNumSolutions()):                        for k in range(self.getNumSolutions()):
2081                           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
2082                        # flux=C-util.reorderComponents(B,[1,0])                        # flux=C-util.reorderComponents(B,[1,0])
2083                          length_of_flux=util.sqrt(flux2)
2084                     else:                     else:
2085                        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)  
2086              elif C.isEmpty():              elif C.isEmpty():
2087                length_of_flux=util.length(B)                length_of_flux=util.length(B)
               #flux=B  
2088              else:              else:
2089                length_of_flux=util.length(C)                length_of_flux=util.length(C)
               #flux=C  
   
             #length_of_flux=util.length(flux)  
2090              flux_max=util.Lsup(length_of_flux)              flux_max=util.Lsup(length_of_flux)
2091              if flux_max>0.:              if flux_max>0.:
2092                 # length_of_A=util.inner(flux,util.tensormutiply(A,flux))                if A.isEmpty():
2093                 length_of_A=util.length(A)                    inv_A=1./self.SMALL_TOLERANCE
2094                 A_max=util.Lsup(length_of_A)                    peclet_number=escript.Scalar(inv_A,length_of_flux.getFunctionSpace())
2095                 if A_max>0:                    xi=self.__xi(self,peclet_number)
2096                      inv_A=1./(length_of_A+A_max*self.__TOL)                else:
2097                 else:                    # length_of_A=util.inner(flux,util.tensormutiply(A,flux))
2098                      inv_A=1./self.__TOL                    length_of_A=util.length(A)
2099                 peclet_number=length_of_flux*h/2*inv_A                    A_max=util.Lsup(length_of_A)
2100                 xi=self.__xi(peclet_number)                    if A_max>0:
2101                 self.__Xi=h*xi/(length_of_flux+flux_max*self.__TOL)                         inv_A=1./(length_of_A+A_max*self.SMALL_TOLERANCE)
2102                 self.trace("preclet number = %e"%util.Lsup(peclet_number))                    else:
2103                           inv_A=1./self.SMALL_TOLERANCE
2104                      peclet_number=length_of_flux*h/2*inv_A
2105                      xi=self.__xi(self,peclet_number)
2106                  self.__Xi=h*xi/(length_of_flux+flux_max*self.SMALL_TOLERANCE)
2107                  self.trace("preclet number = %e"%util.Lsup(peclet_number))
2108                else:
2109                  self.__Xi=escript.Scalar(0.,length_of_flux.getFunctionSpace())
2110        return self.__Xi        return self.__Xi
2111    
2112    
# Line 2101  class AdvectivePDE(LinearPDE): Line 2133  class AdvectivePDE(LinearPDE):
2133              Aout=A              Aout=A
2134           else:           else:
2135              if A.isEmpty():              if A.isEmpty():
2136                 Aout=self.createNewCoefficient("A")                 Aout=self.createCoefficientOfGeneralPDE("A")
2137              else:              else:
2138                 Aout=A[:]                 Aout=A[:]
2139              Xi=self.__getXi()              Xi=self.__getXi()
# Line 2121  class AdvectivePDE(LinearPDE): Line 2153  class AdvectivePDE(LinearPDE):
2153                                 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]
2154                                 # 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)
2155              else:              else:
2156                  for j in range(self.getDim()):                 if not C.isEmpty() and not B.isEmpty():
2157                     for l in range(self.getDim()):                     delta=(C-B)
2158                        if not C.isEmpty() and not B.isEmpty():                     Aout+=util.outer(Xi*delta,delta)
2159                            Aout[j,l]+=Xi*(C[j]-B[j])*(C[l]-B[l])                 elif not B.isEmpty():
2160                        elif C.isEmpty():                     Aout+=util.outer(Xi*B,B)
2161                            Aout[j,l]+=Xi*B[j]*B[l]                 elif not C.isEmpty():
2162                        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)  
2163           return Aout           return Aout
2164       elif name == "B" :       elif name == "B" :
2165             # return self.getCoefficient("B")
2166           B=self.getCoefficient("B")           B=self.getCoefficient("B")
2167           C=self.getCoefficient("C")           C=self.getCoefficient("C")
2168           D=self.getCoefficient("D")           D=self.getCoefficient("D")
# Line 2146  class AdvectivePDE(LinearPDE): Line 2171  class AdvectivePDE(LinearPDE):
2171           else:           else:
2172              Xi=self.__getXi()              Xi=self.__getXi()
2173              if B.isEmpty():              if B.isEmpty():
2174                  Bout=self.createNewCoefficient("B")                  Bout=self.createCoefficientOfGeneralPDE("B")
2175              else:              else:
2176                  Bout=B[:]                  Bout=B[:]
2177              if self.getNumEquations()>1:              if self.getNumEquations()>1:
# Line 2158  class AdvectivePDE(LinearPDE): Line 2183  class AdvectivePDE(LinearPDE):
2183                             Bout[i,j,k]+=tmp*C[p,i,j]                             Bout[i,j,k]+=tmp*C[p,i,j]
2184                             # 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)
2185              else:              else:
2186                 tmp=Xi*D                 Bout+=(Xi*D)*C
                for j in range(self.getDim()): Bout[j]+=tmp*C[j]  
                # Bout=Bout+Xi*D*C  
2187           return Bout           return Bout
2188       elif name == "C" :       elif name == "C" :
2189             # return self.getCoefficient("C")
2190           B=self.getCoefficient("B")           B=self.getCoefficient("B")
2191           C=self.getCoefficient("C")           C=self.getCoefficient("C")
2192           D=self.getCoefficient("D")           D=self.getCoefficient("D")
# Line 2171  class AdvectivePDE(LinearPDE): Line 2195  class AdvectivePDE(LinearPDE):
2195           else:           else:
2196              Xi=self.__getXi()              Xi=self.__getXi()
2197              if C.isEmpty():              if C.isEmpty():
2198                  Cout=self.createNewCoefficient("C")                  Cout=self.createCoefficientOfGeneralPDE("C")
2199              else:              else:
2200                  Cout=C[:]                  Cout=C[:]
2201              if self.getNumEquations()>1:              if self.getNumEquations()>1:
# Line 2183  class AdvectivePDE(LinearPDE): Line 2207  class AdvectivePDE(LinearPDE):
2207                                   Cout[i,k,l]+=tmp*B[p,l,i]                                   Cout[i,k,l]+=tmp*B[p,l,i]
2208                                   # Cout=Cout+Xi*B[p,l,i]*D[p,k]                                   # Cout=Cout+Xi*B[p,l,i]*D[p,k]
2209              else:              else:
2210                 tmp=Xi*D                 Cout+=(Xi*D)*B
                for j in range(self.getDim()): Cout[j]+=tmp*B[j]  
                # Cout=Cout+tmp*D*B  
2211           return Cout           return Cout
2212       elif name == "D" :       elif name == "D" :
2213           return self.getCoefficient("D")           return self.getCoefficient("D")
2214       elif name == "X" :       elif name == "X" :
2215             # return self.getCoefficient("X")
2216           X=self.getCoefficient("X")           X=self.getCoefficient("X")
2217           Y=self.getCoefficient("Y")           Y=self.getCoefficient("Y")
2218           B=self.getCoefficient("B")           B=self.getCoefficient("B")
# Line 2198  class AdvectivePDE(LinearPDE): Line 2221  class AdvectivePDE(LinearPDE):
2221              Xout=X              Xout=X
2222           else:           else:
2223              if X.isEmpty():              if X.isEmpty():
2224                  Xout=self.createNewCoefficient("X")                  Xout=self.createCoefficientOfGeneralPDE("X")
2225              else:              else:
2226                  Xout=X[:]                  Xout=X[:]
2227              Xi=self.__getXi()              Xi=self.__getXi()
# Line 2217  class AdvectivePDE(LinearPDE): Line 2240  class AdvectivePDE(LinearPDE):
2240                               Xout[i,j]+=tmp*C[p,i,j]                               Xout[i,j]+=tmp*C[p,i,j]
2241                               # Xout=X_out+Xi*util.inner(Y,C,offset=1)                               # Xout=X_out+Xi*util.inner(Y,C,offset=1)
2242              else:              else:
2243                   tmp=Xi*Y                if not C.isEmpty() and not B.isEmpty():
2244                   for j in range(self.getDim()):                  Xout+=(Xi*Y)*(C-B)
2245                      if not C.isEmpty() and not B.isEmpty():                elif C.isEmpty():
2246                         Xout[j]+=tmp*(C[j]-B[j])                  Xout-=(Xi*Y)*B
2247                         # Xout=Xout+Xi*Y*(C-B)                else:
2248                      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  
2249           return Xout           return Xout
2250       elif name == "Y" :       elif name == "Y" :
2251           return self.getCoefficient("Y")           return self.getCoefficient("Y")
# Line 2246  class AdvectivePDE(LinearPDE): Line 2264  class AdvectivePDE(LinearPDE):
2264       else:       else:
2265          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
2266    
 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,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}  
   
    @remark: no upwinding is applied yet.  
   
    """  
   
    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(Helmholtz, 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),  
                         "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 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(Helmholtz, 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 escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))*self.getCoefficient("k")  
      elif name == "B" :  
          return escript.Data()  
      elif name == "C" :  
          return escript.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  
   
   
2267  # $Log$  # $Log$
2268  # Revision 1.14  2005/09/22 01:54:57  jgs  # Revision 1.14  2005/09/22 01:54:57  jgs
2269  # 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.345  
changed lines
  Added in v.614

  ViewVC Help
Powered by ViewVC 1.1.26