/[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 596 by gross, Wed Mar 15 22:41:23 2006 UTC revision 969 by ksteube, Tue Feb 13 23:02:23 2007 UTC
# Line 1  Line 1 
1  # $Id$  # $Id$
   
 #  
 #      COPYRIGHT ACcESS 2004 -  All Rights Reserved  
 #  
 #   This software is the property of ACcESS.  No part of this code  
 #   may be copied in any form or by any means without the expressed written  
 #   consent of ACcESS.  Copying, use or modification of this software  
 #   by any unauthorised person is illegal unless that  
 #   person has a software license agreement with ACcESS.  
 #  
2  """  """
3  The module provides an interface to define and solve linear partial  The module provides an interface to define and solve linear partial
4  differential equations (PDEs) within L{escript}. L{linearPDEs} does not provide any  differential equations (PDEs) within L{escript}. L{linearPDEs} does not provide any
# Line 17  the PDE solver library defined through t Line 7  the PDE solver library defined through t
7  The general interface is provided through the L{LinearPDE} class. The  The general interface is provided through the L{LinearPDE} class. The
8  L{AdvectivePDE} which is derived from the L{LinearPDE} class  L{AdvectivePDE} which is derived from the L{LinearPDE} class
9  provides an interface to PDE dominated by its advective terms. The L{Poisson},  provides an interface to PDE dominated by its advective terms. The L{Poisson},
10  L{Helmholtz}, L{LameEquation}, L{AdvectionDiffusion}  L{Helmholtz}, L{LameEquation}, L{AdvectivePDE}
11  classs which are also derived form the L{LinearPDE} class should be used  classs which are also derived form the L{LinearPDE} class should be used
12  to define of solve these sepecial PDEs.  to define of solve these sepecial PDEs.
13    
14  @var __author__: name of author  @var __author__: name of author
15  @var __licence__: licence agreement  @var __copyright__: copyrights
16    @var __license__: licence agreement
17  @var __url__: url entry point on documentation  @var __url__: url entry point on documentation
18  @var __version__: version  @var __version__: version
19  @var __date__: date of the version  @var __date__: date of the version
# Line 33  import util Line 24  import util
24  import numarray  import numarray
25    
26  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
27  __licence__="contact: esys@access.uq.edu.au"  __copyright__="""  Copyright (c) 2006 by ACcESS MNRF
28  __url__="http://www.iservo.edu.au/esys/escript"                      http://www.access.edu.au
29                    Primary Business: Queensland, Australia"""
30    __license__="""Licensed under the Open Software License version 3.0
31                 http://www.opensource.org/licenses/osl-3.0.php"""
32    __url__="http://www.iservo.edu.au/esys"
33  __version__="$Revision$"  __version__="$Revision$"
34  __date__="$Date$"  __date__="$Date$"
35    
# Line 120  class PDECoefficient(object): Line 115  class PDECoefficient(object):
115         @param domain: domain on which the PDE uses the coefficient         @param domain: domain on which the PDE uses the coefficient
116         @type domain: L{Domain<escript.Domain>}         @type domain: L{Domain<escript.Domain>}
117         @param reducedEquationOrder: True to indicate that reduced order is used to represent the equation         @param reducedEquationOrder: True to indicate that reduced order is used to represent the equation
118         @type domain: C{bool}         @type reducedEquationOrder: C{bool}
119         @param reducedSolutionOrder: True to indicate that reduced order is used to represent the solution         @param reducedSolutionOrder: True to indicate that reduced order is used to represent the solution
120         @type domain: C{bool}         @type reducedSolutionOrder: C{bool}
121         @return:  L{FunctionSpace<escript.FunctionSpace>} of the coefficient         @return:  L{FunctionSpace<escript.FunctionSpace>} of the coefficient
122         @rtype:  L{FunctionSpace<escript.FunctionSpace>}         @rtype:  L{FunctionSpace<escript.FunctionSpace>}
123         """         """
# Line 160  class PDECoefficient(object): Line 155  class PDECoefficient(object):
155         @param numSolutions: number of components of the PDE solution         @param numSolutions: number of components of the PDE solution
156         @type numSolutions: C{int}         @type numSolutions: C{int}
157         @param reducedEquationOrder: True to indicate that reduced order is used to represent the equation         @param reducedEquationOrder: True to indicate that reduced order is used to represent the equation
158         @type domain: C{bool}         @type reducedEquationOrder: C{bool}
159         @param reducedSolutionOrder: True to indicate that reduced order is used to represent the solution         @param reducedSolutionOrder: True to indicate that reduced order is used to represent the solution
160         @type domain: C{bool}         @type reducedSolutionOrder: C{bool}
161         @param newValue: number of components of the PDE solution         @param newValue: number of components of the PDE solution
162         @type newValue: any object that can be converted into a L{Data<escript.Data>} object with the appropriate shape and L{FunctionSpace<escript.FunctionSpace>}         @type newValue: any object that can be converted into a L{Data<escript.Data>} object with the appropriate shape and L{FunctionSpace<escript.FunctionSpace>}
163         @raise IllegalCoefficientValue: if the shape of the assigned value does not match the shape of the coefficient         @raise IllegalCoefficientValue: if the shape of the assigned value does not match the shape of the coefficient
# Line 418  class LinearPDE(object): Line 413  class LinearPDE(object):
413     @cvar SCSL: SGI SCSL solver library     @cvar SCSL: SGI SCSL solver library
414     @cvar MKL: Intel's MKL solver library     @cvar MKL: Intel's MKL solver library
415     @cvar UMFPACK: the UMFPACK library     @cvar UMFPACK: the UMFPACK library
416       @cvar TRILINOS: the TRILINOS parallel solver class library from Sandia Natl Labs
417     @cvar ITERATIVE: The default iterative solver     @cvar ITERATIVE: The default iterative solver
418     @cvar AMG: algebraic multi grid     @cvar AMG: algebraic multi grid
419     @cvar RILU: recursive ILU     @cvar RILU: recursive ILU
# Line 447  class LinearPDE(object): Line 443  class LinearPDE(object):
443     PASO= 21     PASO= 21
444     AMG= 22     AMG= 22
445     RILU = 23     RILU = 23
446       TRILINOS = 24
447    
448     SMALL_TOLERANCE=1.e-13     SMALL_TOLERANCE=1.e-13
449     __PACKAGE_KEY="package"     __PACKAGE_KEY="package"
# Line 669  class LinearPDE(object): Line 666  class LinearPDE(object):
666       @rtype: L{Data<escript.Data>}       @rtype: L{Data<escript.Data>}
667       """       """
668       if u==None:       if u==None:
669            return self.getOperator()*self.getSolution()          return self.getOperator()*self.getSolution()
670       else:       else:
671          self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())          return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())
672    
673     def getResidual(self,u=None):     def getResidual(self,u=None):
674       """       """
# Line 873  class LinearPDE(object): Line 870  class LinearPDE(object):
870         elif p==self.MKL: package= "MKL"         elif p==self.MKL: package= "MKL"
871         elif p==self.SCSL: package= "SCSL"         elif p==self.SCSL: package= "SCSL"
872         elif p==self.UMFPACK: package= "UMFPACK"         elif p==self.UMFPACK: package= "UMFPACK"
873           elif p==self.TRILINOS: package= "TRILINOS"
874         else : method="unknown"         else : method="unknown"
875         return "%s solver of %s package"%(method,package)         return "%s solver of %s package"%(method,package)
876    
# Line 890  class LinearPDE(object): Line 888  class LinearPDE(object):
888         """         """
889         sets a new solver package         sets a new solver package
890    
891         @param solver: sets a new solver method.         @param package: sets a new solver method.
892         @type solver: one of L{DEFAULT}, L{PASO} L{SCSL}, L{MKL}, L{UMLPACK}         @type package: one of L{DEFAULT}, L{PASO} L{SCSL}, L{MKL}, L{UMFPACK}, L{TRILINOS}
893         """         """
894         if package==None: package=self.DEFAULT         if package==None: package=self.DEFAULT
895         if not package==self.getSolverPackage():         if not package==self.getSolverPackage():
896             self.__solver_method=solver             self.__solver_package=package
897             self.__checkMatrixType()             self.__checkMatrixType()
898             self.trace("New solver is %s"%self.getSolverMethodName())             self.trace("New solver is %s"%self.getSolverMethodName())
899    
# Line 928  class LinearPDE(object): Line 926  class LinearPDE(object):
926         @param tol: new tolerance for the solver. If the tol is lower then the current tolerence         @param tol: new tolerance for the solver. If the tol is lower then the current tolerence
927                     the system will be resolved.                     the system will be resolved.
928         @type tol: positive C{float}         @type tol: positive C{float}
929         @raise ValueException: if tolerance is not positive.         @raise ValueError: if tolerance is not positive.
930         """         """
931         if not tol>0:         if not tol>0:
932             raise ValueException,"Tolerance as to be positive"             raise ValueError,"Tolerance as to be positive"
933         if tol<self.getTolerance(): self.__invalidateSolution()         if tol<self.getTolerance(): self.__invalidateSolution()
934         self.trace("New tolerance %e"%tol)         self.trace("New tolerance %e"%tol)
935         self.__tolerance=tol         self.__tolerance=tol
# Line 1518  class LinearPDE(object): Line 1516  class LinearPDE(object):
1516           r=self.getCoefficientOfGeneralPDE("r")           r=self.getCoefficientOfGeneralPDE("r")
1517           homogeneous_constraint=True           homogeneous_constraint=True
1518           if not q.isEmpty() and not r.isEmpty():           if not q.isEmpty() and not r.isEmpty():
1519               if util.Lsup(q*r)>=1.e-13*util.Lsup(r):               if util.Lsup(q*r)>0.:
1520                 self.trace("Inhomogeneous constraint detected.")                 self.trace("Inhomogeneous constraint detected.")
1521                 self.__invalidateSystem()                 self.__invalidateSystem()
1522    
# Line 1543  class LinearPDE(object): Line 1541  class LinearPDE(object):
1541                   D=self.getCoefficientOfGeneralPDE("D")                   D=self.getCoefficientOfGeneralPDE("D")
1542                   if not D.isEmpty():                   if not D.isEmpty():
1543                       if self.getNumSolutions()>1:                       if self.getNumSolutions()>1:
1544                          D_times_e=util.matrixmult(D,numarray.ones((self.getNumSolutions(),)))                          #D_times_e=util.matrixmult(D,numarray.ones((self.getNumSolutions(),)))
1545                            D_times_e=util.matrix_mult(D,numarray.ones((self.getNumSolutions(),)))
1546                       else:                       else:
1547                          D_times_e=D                          D_times_e=D
1548                   else:                   else:
# Line 1551  class LinearPDE(object): Line 1550  class LinearPDE(object):
1550                   d=self.getCoefficientOfGeneralPDE("d")                   d=self.getCoefficientOfGeneralPDE("d")
1551                   if not d.isEmpty():                   if not d.isEmpty():
1552                       if self.getNumSolutions()>1:                       if self.getNumSolutions()>1:
1553                          d_times_e=util.matrixmult(d,numarray.ones((self.getNumSolutions(),)))                          #d_times_e=util.matrixmult(d,numarray.ones((self.getNumSolutions(),)))
1554                            d_times_e=util.matrix_mult(d,numarray.ones((self.getNumSolutions(),)))
1555                       else:                       else:
1556                          d_times_e=d                          d_times_e=d
1557                   else:                   else:
# Line 1812  class LameEquation(LinearPDE): Line 1812  class LameEquation(LinearPDE):
1812     """     """
1813     Class to define a Lame equation problem:     Class to define a Lame equation problem:
1814    
1815     M{-grad(S{mu}*(grad(u[i])[j]+grad(u[j])[i]))[j] - grad(S{lambda}*grad(u[j])[i])[j] = F_i -grad(S{sigma}[i,j])[j] }     M{-grad(S{mu}*(grad(u[i])[j]+grad(u[j])[i]))[j] - grad(S{lambda}*grad(u[k])[k])[j] = F_i -grad(S{sigma}[ij])[j] }
1816    
1817     with natural boundary conditons:     with natural boundary conditons:
1818    
1819     M{n[j]*(S{mu}*(grad(u[i])[j]+grad(u[j])[i]) - S{lambda}*grad(u[j])[i]) = f_i -n[j]*S{sigma}[i,j] }     M{n[j]*(S{mu}*(grad(u[i])[j]+grad(u[j])[i]) + S{lambda}*grad(u[k])[k]) = f_i +n[j]*S{sigma}[ij] }
1820    
1821     and constraints:     and constraints:
1822    

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

  ViewVC Help
Powered by ViewVC 1.1.26