/[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 328 by gross, Wed Dec 7 04:41:53 2005 UTC revision 637 by gross, Thu Mar 23 10:55:31 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{Helmholtz}, L{LameEquation}, L{AdvectionDiffusion}
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 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"
450     __TOLERANCE_KEY="tolerance"     __TOLERANCE_KEY="tolerance"
451       __PRECONDITIONER_KEY="preconditioner"
452    
453    
454     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 498  class LinearPDE(object):
498       self.__tolerance=1.e-8       self.__tolerance=1.e-8
499       self.__solver_method=self.DEFAULT       self.__solver_method=self.DEFAULT
500       self.__solver_package=self.DEFAULT       self.__solver_package=self.DEFAULT
501         self.__preconditioner=self.DEFAULT
502       self.__matrix_type=self.__domain.getSystemMatrixTypeId(self.DEFAULT,self.DEFAULT,False)       self.__matrix_type=self.__domain.getSystemMatrixTypeId(self.DEFAULT,self.DEFAULT,False)
503       self.__sym=False       self.__sym=False
504    
# Line 697  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 721  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 737  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 745  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 753  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 772  class LinearPDE(object): Line 773  class LinearPDE(object):
773         @type verbose: C{bool}         @type verbose: C{bool}
774         @keyword reordering: reordering scheme to be used during elimination. Allowed values are         @keyword reordering: reordering scheme to be used during elimination. Allowed values are
775                              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}  
776         @keyword iter_max: maximum number of iteration steps allowed.         @keyword iter_max: maximum number of iteration steps allowed.
777         @keyword drop_tolerance: threshold for drupping in L{ILUT}         @keyword drop_tolerance: threshold for drupping in L{ILUT}
778         @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 785  class LinearPDE(object):
785               self.__solution=self.copyConstraint(f*mat)               self.__solution=self.copyConstraint(f*mat)
786            else:            else:
787               options[self.__TOLERANCE_KEY]=self.getTolerance()               options[self.__TOLERANCE_KEY]=self.getTolerance()
788               options[self.__METHOD_KEY]=self.getSolverMethod()               options[self.__METHOD_KEY]=self.getSolverMethod()[0]
789                 options[self.__PRECONDITIONER_KEY]=self.getSolverMethod()[1]
790               options[self.__PACKAGE_KEY]=self.getSolverPackage()               options[self.__PACKAGE_KEY]=self.getSolverPackage()
791               options[self.__SYMMETRY_KEY]=self.isSymmetric()               options[self.__SYMMETRY_KEY]=self.isSymmetric()
792               self.trace("PDE is resolved.")               self.trace("PDE is resolved.")
# Line 815  class LinearPDE(object): Line 815  class LinearPDE(object):
815     # =============================================================================     # =============================================================================
816     #   solver settings:     #   solver settings:
817     # =============================================================================     # =============================================================================
818     def setSolverMethod(self,solver=None):     def setSolverMethod(self,solver=None,preconditioner=None):
819         """         """
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.
825           @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 not solver==self.getSolverMethod():         if preconditioner==None: preconditioner=self.DEFAULT
829           if not (solver,preconditioner)==self.getSolverMethod():
830             self.__solver_method=solver             self.__solver_method=solver
831               self.__preconditioner=preconditioner
832             self.__checkMatrixType()             self.__checkMatrixType()
833             self.trace("New solver is %s"%self.getSolverMethodName())             self.trace("New solver is %s"%self.getSolverMethodName())
834    
# Line 838  class LinearPDE(object): Line 842  class LinearPDE(object):
842    
843         m=self.getSolverMethod()         m=self.getSolverMethod()
844         p=self.getSolverPackage()         p=self.getSolverPackage()
845         if m==self.DEFAULT: method="DEFAULT"         method=""
846         elif m==self.DIRECT: method= "DIRECT"         if m[0]==self.DEFAULT: method="DEFAULT"
847         elif m==self.ITERATIVE: method= "ITERATIVE"         elif m[0]==self.DIRECT: method= "DIRECT"
848         elif m==self.CHOLEVSKY: method= "CHOLEVSKY"         elif m[0]==self.ITERATIVE: method= "ITERATIVE"
849         elif m==self.PCG: method= "PCG"         elif m[0]==self.CHOLEVSKY: method= "CHOLEVSKY"
850         elif m==self.CR: method= "CR"         elif m[0]==self.PCG: method= "PCG"
851         elif m==self.CGS: method= "CGS"         elif m[0]==self.CR: method= "CR"
852         elif m==self.BICGSTAB: method= "BICGSTAB"         elif m[0]==self.CGS: method= "CGS"
853         elif m==self.SSOR: method= "SSOR"         elif m[0]==self.BICGSTAB: method= "BICGSTAB"
854         elif m==self.GMRES: method= "GMRES"         elif m[0]==self.SSOR: method= "SSOR"
855         elif m==self.PRES20: method= "PRES20"         elif m[0]==self.GMRES: method= "GMRES"
856         elif m==self.LUMPING: method= "LUMPING"         elif m[0]==self.PRES20: method= "PRES20"
857         else : method="unknown"         elif m[0]==self.LUMPING: method= "LUMPING"
858           elif m[0]==self.AMG: method= "AMG"
859           if m[1]==self.DEFAULT: method+="+DEFAULT"
860           elif m[1]==self.JACOBI: method+= "+JACOBI"
861           elif m[1]==self.ILU0: method+= "+ILU0"
862           elif m[1]==self.ILUT: method+= "+ILUT"
863           elif m[1]==self.SSOR: method+= "+SSOR"
864           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 867  class LinearPDE(object): Line 879  class LinearPDE(object):
879         @return: the solver method currently be used.         @return: the solver method currently be used.
880         @rtype: C{int}         @rtype: C{int}
881         """         """
882         return self.__solver_method         return self.__solver_method,self.__preconditioner
883    
884     def setSolverPackage(self,package=None):     def setSolverPackage(self,package=None):
885         """         """
# Line 898  class LinearPDE(object): Line 910  class LinearPDE(object):
910        @return: True is lumping is currently used a solver method.        @return: True is lumping is currently used a solver method.
911        @rtype: C{bool}        @rtype: C{bool}
912        """        """
913        return self.getSolverMethod()==self.LUMPING        return self.getSolverMethod()[0]==self.LUMPING
914    
915     def setTolerance(self,tol=1.e-8):     def setTolerance(self,tol=1.e-8):
916         """         """
# Line 1091  class LinearPDE(object): Line 1103  class LinearPDE(object):
1103       """       """
1104       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.
1105       """       """
1106       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())
1107       if not new_matrix_type==self.__matrix_type:       if not new_matrix_type==self.__matrix_type:
1108           self.trace("Matrix type is now %d."%new_matrix_type)           self.trace("Matrix type is now %d."%new_matrix_type)
1109           self.__matrix_type=new_matrix_type           self.__matrix_type=new_matrix_type
# Line 1515  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"), \                       else:
1543                             escript.Data(),\                          D_times_e=D
1544                             self.getCoefficientOfGeneralPDE("d_contact"), \                   else:
1545                             escript.Data())                      D_times_e=escript.Data()
1546                   self.__operator=1./(mat*escript.Data(1,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True))                   d=self.getCoefficientOfGeneralPDE("d")
1547                   del mat                   if not d.isEmpty():
1548                         if self.getNumSolutions()>1:
1549                            d_times_e=util.matrixmult(d,numarray.ones((self.getNumSolutions(),)))
1550                         else:
1551                            d_times_e=d
1552                     else:
1553                        d_times_e=escript.Data()
1554                     d_contact=self.getCoefficientOfGeneralPDE("d_contact")
1555                     if not d_contact.isEmpty():
1556                         if self.getNumSolutions()>1:
1557                            d_contact_times_e=util.matrixmult(d_contact,numarray.ones((self.getNumSolutions(),)))
1558                         else:
1559                            d_contact_times_e=d_contact
1560                     else:
1561                        d_contact_times_e=escript.Data()
1562        
1563                     self.__operator=self.__getNewRightHandSide()
1564                     self.getDomain().addPDEToRHS(self.__operator, \
1565                                                  escript.Data(), \
1566                                                  D_times_e, \
1567                                                  d_times_e,\
1568                                                  d_contact_times_e)
1569                     self.__operator=1./self.__operator
1570                   self.trace("New lumped operator has been built.")                   self.trace("New lumped operator has been built.")
1571                   self.__operator_is_Valid=True                   self.__operator_is_Valid=True
1572                if not self.__righthandside_isValid:                if not self.__righthandside_isValid:
# Line 1797  class LameEquation(LinearPDE): Line 1831  class LameEquation(LinearPDE):
1831                            "q"            : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}                            "q"            : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}
1832        self.setSymmetryOn()        self.setSymmetryOn()
1833    
1834     def setValue(self,**coefficients):     def setValues(self,**coefficients):
1835       """       """
1836       sets new values to coefficients       sets new values to coefficients
1837    
# Line 1820  class LameEquation(LinearPDE): Line 1854  class LameEquation(LinearPDE):
1854                 depending of reduced order is used for the representation of the equation.                 depending of reduced order is used for the representation of the equation.
1855       @raise IllegalCoefficient: if an unknown coefficient keyword is used.       @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1856       """       """
1857       super(LameEquation, self).setValue(**coefficients)       super(LameEquation, self).setValues(**coefficients)
1858    
1859     def getCoefficientOfGeneralPDE(self,name):     def getCoefficientOfGeneralPDE(self,name):
1860       """       """
# Line 1934  class AdvectivePDE(LinearPDE): Line 1968  class AdvectivePDE(LinearPDE):
1968           self.__xi=xi           self.__xi=xi
1969        self.__Xi=escript.Data()        self.__Xi=escript.Data()
1970    
1971     def setValue(**coefficients):     def setValue(self,**coefficients):
1972        """        """
1973        sets new values to coefficients        sets new values to coefficients
1974    
# Line 2018  class AdvectivePDE(LinearPDE): Line 2052  class AdvectivePDE(LinearPDE):
2052       """       """
2053       return escript.Scalar(0.5,P.getFunctionSpace())       return escript.Scalar(0.5,P.getFunctionSpace())
2054    
    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.  
   
2055     def __getXi(self):     def __getXi(self):
2056        if self.__Xi.isEmpty():        if self.__Xi.isEmpty():
2057           B=self.getCoefficient("B")           B=self.getCoefficient("B")
# Line 2034  class AdvectivePDE(LinearPDE): Line 2061  class AdvectivePDE(LinearPDE):
2061           self.__Xi=escript.Scalar(0.,self.getFunctionSpaceForCoefficient("A"))           self.__Xi=escript.Scalar(0.,self.getFunctionSpaceForCoefficient("A"))
2062           if not C.isEmpty() or not B.isEmpty():           if not C.isEmpty() or not B.isEmpty():
2063              if not C.isEmpty() and not B.isEmpty():              if not C.isEmpty() and not B.isEmpty():
                 flux2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))  
2064                  if self.getNumEquations()>1:                  if self.getNumEquations()>1:
2065                     if self.getNumSolutions()>1:                     if self.getNumSolutions()>1:
2066                          flux2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))
2067                        for i in range(self.getNumEquations()):                        for i in range(self.getNumEquations()):
2068                           for k in range(self.getNumSolutions()):                           for k in range(self.getNumSolutions()):
2069                              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
2070                          length_of_flux=util.sqrt(flux2)
2071                        # flux=C-util.reorderComponents(B,[0,2,1])                        # flux=C-util.reorderComponents(B,[0,2,1])
2072                     else:                     else:
2073                          flux2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))
2074                        for i in range(self.getNumEquations()):                        for i in range(self.getNumEquations()):
2075                           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
2076                          length_of_flux=util.sqrt(flux2)
2077                        # flux=C-B                        # flux=C-B
2078                  else:                  else:
2079                     if self.getNumSolutions()>1:                     if self.getNumSolutions()>1:
2080                          flux2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))
2081                        for k in range(self.getNumSolutions()):                        for k in range(self.getNumSolutions()):
2082                           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
2083                        # flux=C-util.reorderComponents(B,[1,0])                        # flux=C-util.reorderComponents(B,[1,0])
2084                          length_of_flux=util.sqrt(flux2)
2085                     else:                     else:
2086                        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)  
2087              elif C.isEmpty():              elif C.isEmpty():
2088                length_of_flux=util.length(B)                length_of_flux=util.length(B)
               #flux=B  
2089              else:              else:
2090                length_of_flux=util.length(C)                length_of_flux=util.length(C)
               #flux=C  
   
             #length_of_flux=util.length(flux)  
2091              flux_max=util.Lsup(length_of_flux)              flux_max=util.Lsup(length_of_flux)
2092              if flux_max>0.:              if flux_max>0.:
2093                 # length_of_A=util.inner(flux,util.tensormutiply(A,flux))                if A.isEmpty():
2094                 length_of_A=util.length(A)                    inv_A=1./self.SMALL_TOLERANCE
2095                 A_max=util.Lsup(length_of_A)                    peclet_number=escript.Scalar(inv_A,length_of_flux.getFunctionSpace())
2096                 if A_max>0:                    xi=self.__xi(self,peclet_number)
2097                      inv_A=1./(length_of_A+A_max*self.__TOL)                else:
2098                 else:                    # length_of_A=util.inner(flux,util.tensormutiply(A,flux))
2099                      inv_A=1./self.__TOL                    length_of_A=util.length(A)
2100                 peclet_number=length_of_flux*h/2*inv_A                    A_max=util.Lsup(length_of_A)
2101                 xi=self.__xi(peclet_number)                    if A_max>0:
2102                 self.__Xi=h*xi/(length_of_flux+flux_max*self.__TOL)                         inv_A=1./(length_of_A+A_max*self.SMALL_TOLERANCE)
2103                 self.trace("preclet number = %e"%util.Lsup(peclet_number))                    else:
2104                           inv_A=1./self.SMALL_TOLERANCE
2105                      peclet_number=length_of_flux*h/2*inv_A
2106                      xi=self.__xi(self,peclet_number)
2107                  self.__Xi=h*xi/(length_of_flux+flux_max*self.SMALL_TOLERANCE)
2108                  self.trace("preclet number = %e"%util.Lsup(peclet_number))
2109                else:
2110                  self.__Xi=escript.Scalar(0.,length_of_flux.getFunctionSpace())
2111        return self.__Xi        return self.__Xi
2112    
2113    
# Line 2101  class AdvectivePDE(LinearPDE): Line 2134  class AdvectivePDE(LinearPDE):
2134              Aout=A              Aout=A
2135           else:           else:
2136              if A.isEmpty():              if A.isEmpty():
2137                 Aout=self.createNewCoefficient("A")                 Aout=self.createCoefficientOfGeneralPDE("A")
2138              else:              else:
2139                 Aout=A[:]                 Aout=A[:]
2140              Xi=self.__getXi()              Xi=self.__getXi()
# Line 2121  class AdvectivePDE(LinearPDE): Line 2154  class AdvectivePDE(LinearPDE):
2154                                 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]
2155                                 # 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)
2156              else:              else:
2157                  for j in range(self.getDim()):                 if not C.isEmpty() and not B.isEmpty():
2158                     for l in range(self.getDim()):                     delta=(C-B)
2159                        if not C.isEmpty() and not B.isEmpty():                     Aout+=util.outer(Xi*delta,delta)
2160                            Aout[j,l]+=Xi*(C[j]-B[j])*(C[l]-B[l])                 elif not B.isEmpty():
2161                        elif C.isEmpty():                     Aout+=util.outer(Xi*B,B)
2162                            Aout[j,l]+=Xi*B[j]*B[l]                 elif not C.isEmpty():
2163                        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)  
2164           return Aout           return Aout
2165       elif name == "B" :       elif name == "B" :
2166             # return self.getCoefficient("B")
2167           B=self.getCoefficient("B")           B=self.getCoefficient("B")
2168           C=self.getCoefficient("C")           C=self.getCoefficient("C")
2169           D=self.getCoefficient("D")           D=self.getCoefficient("D")
# Line 2146  class AdvectivePDE(LinearPDE): Line 2172  class AdvectivePDE(LinearPDE):
2172           else:           else:
2173              Xi=self.__getXi()              Xi=self.__getXi()
2174              if B.isEmpty():              if B.isEmpty():
2175                  Bout=self.createNewCoefficient("B")                  Bout=self.createCoefficientOfGeneralPDE("B")
2176              else:              else:
2177                  Bout=B[:]                  Bout=B[:]
2178              if self.getNumEquations()>1:              if self.getNumEquations()>1:
# Line 2158  class AdvectivePDE(LinearPDE): Line 2184  class AdvectivePDE(LinearPDE):
2184                             Bout[i,j,k]+=tmp*C[p,i,j]                             Bout[i,j,k]+=tmp*C[p,i,j]
2185                             # 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)
2186              else:              else:
2187                 tmp=Xi*D                 Bout+=(Xi*D)*C
                for j in range(self.getDim()): Bout[j]+=tmp*C[j]  
                # Bout=Bout+Xi*D*C  
2188           return Bout           return Bout
2189       elif name == "C" :       elif name == "C" :
2190             # return self.getCoefficient("C")
2191           B=self.getCoefficient("B")           B=self.getCoefficient("B")
2192           C=self.getCoefficient("C")           C=self.getCoefficient("C")
2193           D=self.getCoefficient("D")           D=self.getCoefficient("D")
# Line 2171  class AdvectivePDE(LinearPDE): Line 2196  class AdvectivePDE(LinearPDE):
2196           else:           else:
2197              Xi=self.__getXi()              Xi=self.__getXi()
2198              if C.isEmpty():              if C.isEmpty():
2199                  Cout=self.createNewCoefficient("C")                  Cout=self.createCoefficientOfGeneralPDE("C")
2200              else:              else:
2201                  Cout=C[:]                  Cout=C[:]
2202              if self.getNumEquations()>1:              if self.getNumEquations()>1:
# Line 2183  class AdvectivePDE(LinearPDE): Line 2208  class AdvectivePDE(LinearPDE):
2208                                   Cout[i,k,l]+=tmp*B[p,l,i]                                   Cout[i,k,l]+=tmp*B[p,l,i]
2209                                   # Cout=Cout+Xi*B[p,l,i]*D[p,k]                                   # Cout=Cout+Xi*B[p,l,i]*D[p,k]
2210              else:              else:
2211                 tmp=Xi*D                 Cout+=(Xi*D)*B
                for j in range(self.getDim()): Cout[j]+=tmp*B[j]  
                # Cout=Cout+tmp*D*B  
2212           return Cout           return Cout
2213       elif name == "D" :       elif name == "D" :
2214           return self.getCoefficient("D")           return self.getCoefficient("D")
2215       elif name == "X" :       elif name == "X" :
2216             # return self.getCoefficient("X")
2217           X=self.getCoefficient("X")           X=self.getCoefficient("X")
2218           Y=self.getCoefficient("Y")           Y=self.getCoefficient("Y")
2219           B=self.getCoefficient("B")           B=self.getCoefficient("B")
# Line 2198  class AdvectivePDE(LinearPDE): Line 2222  class AdvectivePDE(LinearPDE):
2222              Xout=X              Xout=X
2223           else:           else:
2224              if X.isEmpty():              if X.isEmpty():
2225                  Xout=self.createNewCoefficient("X")                  Xout=self.createCoefficientOfGeneralPDE("X")
2226              else:              else:
2227                  Xout=X[:]                  Xout=X[:]
2228              Xi=self.__getXi()              Xi=self.__getXi()
# Line 2217  class AdvectivePDE(LinearPDE): Line 2241  class AdvectivePDE(LinearPDE):
2241                               Xout[i,j]+=tmp*C[p,i,j]                               Xout[i,j]+=tmp*C[p,i,j]
2242                               # Xout=X_out+Xi*util.inner(Y,C,offset=1)                               # Xout=X_out+Xi*util.inner(Y,C,offset=1)
2243              else:              else:
2244                   tmp=Xi*Y                if not C.isEmpty() and not B.isEmpty():
2245                   for j in range(self.getDim()):                  Xout+=(Xi*Y)*(C-B)
2246                      if not C.isEmpty() and not B.isEmpty():                elif C.isEmpty():
2247                         Xout[j]+=tmp*(C[j]-B[j])                  Xout-=(Xi*Y)*B
2248                         # Xout=Xout+Xi*Y*(C-B)                else:
2249                      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  
2250           return Xout           return Xout
2251       elif name == "Y" :       elif name == "Y" :
2252           return self.getCoefficient("Y")           return self.getCoefficient("Y")
# Line 2246  class AdvectivePDE(LinearPDE): Line 2265  class AdvectivePDE(LinearPDE):
2265       else:       else:
2266          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
2267    
   
2268  # $Log$  # $Log$
2269  # Revision 1.14  2005/09/22 01:54:57  jgs  # Revision 1.14  2005/09/22 01:54:57  jgs
2270  # 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.328  
changed lines
  Added in v.637
  ViewVC Help
Powered by ViewVC 1.1.26