# Diff of /trunk/escript/py_src/linearPDEs.py

revision 1137 by gross, Thu May 10 08:11:31 2007 UTC revision 1388 by trankine, Fri Jan 11 07:45:58 2008 UTC
# Line 1  Line 1
1    #
2  # \$Id\$  # \$Id\$
3    #
4    #######################################################
5    #
6    #           Copyright 2003-2007 by ACceSS MNRF
7    #       Copyright 2007 by University of Queensland
8    #
9    #                http://esscc.uq.edu.au
10    #        Primary Business: Queensland, Australia
13    #
14    #######################################################
15    #
16
17  """  """
18  The module provides an interface to define and solve linear partial  The module provides an interface to define and solve linear partial
19  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 112  class PDECoefficient(object): Line 127  class PDECoefficient(object):
127         @param reduced: indicates if reduced         @param reduced: indicates if reduced
128         @type reduced: C{bool}         @type reduced: C{bool}
129         """         """

130         super(PDECoefficient, self).__init__()         super(PDECoefficient, self).__init__()
131         self.what=where         self.what=where
132         self.pattern=pattern         self.pattern=pattern
# Line 365  class LinearPDE(object): Line 379  class LinearPDE(object):
379
380     The PDE is symmetrical if     The PDE is symmetrical if
381
382     M{A[i,j]=A[j,i]}  and M{B[j]=C[j]} and M{A_reduced[i,j]=A_reduced[j,i]}  and M{B_reduced[j]=C_reduced[j]     M{A[i,j]=A[j,i]}  and M{B[j]=C[j]} and M{A_reduced[i,j]=A_reduced[j,i]}  and M{B_reduced[j]=C_reduced[j]}
383
384     For a system of PDEs and a solution with several components the PDE has the form     For a system of PDEs and a solution with several components the PDE has the form
385
# Line 416  class LinearPDE(object): Line 430  class LinearPDE(object):
430     of the solution at side 1 and at side 0, denotes the jump of M{u} across discontinuity along the normal calcualted by     of the solution at side 1 and at side 0, denotes the jump of M{u} across discontinuity along the normal calcualted by
431     L{jump<util.jump>}.     L{jump<util.jump>}.
432     The coefficient M{d_contact} is a rank two and M{y_contact} is a rank one both in the L{FunctionOnContactZero<escript.FunctionOnContactZero>} or L{FunctionOnContactOne<escript.FunctionOnContactOne>}.     The coefficient M{d_contact} is a rank two and M{y_contact} is a rank one both in the L{FunctionOnContactZero<escript.FunctionOnContactZero>} or L{FunctionOnContactOne<escript.FunctionOnContactOne>}.
433      The coefficient M{d_contact_reduced} is a rank two and M{y_contact_reduced} is a rank one both in the L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}.     The coefficient M{d_contact_reduced} is a rank two and M{y_contact_reduced} is a rank one both in the L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}.
434     In case of a single PDE and a single component solution the contact condition takes the form     In case of a single PDE and a single component solution the contact condition takes the form
435
436     M{n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)}     M{n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)}
# Line 444  class LinearPDE(object): Line 458  class LinearPDE(object):
458     @cvar SCSL: SGI SCSL solver library     @cvar SCSL: SGI SCSL solver library
459     @cvar MKL: Intel's MKL solver library     @cvar MKL: Intel's MKL solver library
460     @cvar UMFPACK: the UMFPACK library     @cvar UMFPACK: the UMFPACK library
461       @cvar TRILINOS: the TRILINOS parallel solver class library from Sandia Natl Labs
462     @cvar ITERATIVE: The default iterative solver     @cvar ITERATIVE: The default iterative solver
463     @cvar AMG: algebraic multi grid     @cvar AMG: algebraic multi grid
464     @cvar RILU: recursive ILU     @cvar RILU: recursive ILU
# Line 473  class LinearPDE(object): Line 488  class LinearPDE(object):
488     PASO= 21     PASO= 21
489     AMG= 22     AMG= 22
490     RILU = 23     RILU = 23
491       TRILINOS = 24
492
493     SMALL_TOLERANCE=1.e-13     SMALL_TOLERANCE=1.e-13
494     __PACKAGE_KEY="package"     __PACKAGE_KEY="package"
# Line 935  class LinearPDE(object): Line 951  class LinearPDE(object):
951         @param preconditioner: sets a new solver method.         @param preconditioner: sets a new solver method.
952         @type preconditioner: one of L{DEFAULT}, L{JACOBI} L{ILU0}, L{ILUT},L{SSOR}, L{RILU}         @type preconditioner: one of L{DEFAULT}, L{JACOBI} L{ILU0}, L{ILUT},L{SSOR}, L{RILU}
953         """         """
954         if solver==None: solve=self.DEFAULT         if solver==None: solver=self.__solver_method
955           if preconditioner==None: preconditioner=self.__preconditioner
956           if solver==None: solver=self.DEFAULT
957         if preconditioner==None: preconditioner=self.DEFAULT         if preconditioner==None: preconditioner=self.DEFAULT
958         if not (solver,preconditioner)==self.getSolverMethod():         if not (solver,preconditioner)==self.getSolverMethod():
959             self.__solver_method=solver             self.__solver_method=solver
# Line 979  class LinearPDE(object): Line 997  class LinearPDE(object):
997         elif p==self.MKL: package= "MKL"         elif p==self.MKL: package= "MKL"
998         elif p==self.SCSL: package= "SCSL"         elif p==self.SCSL: package= "SCSL"
999         elif p==self.UMFPACK: package= "UMFPACK"         elif p==self.UMFPACK: package= "UMFPACK"
1000           elif p==self.TRILINOS: package= "TRILINOS"
1001         else : method="unknown"         else : method="unknown"
1002         return "%s solver of %s package"%(method,package)         return "%s solver of %s package"%(method,package)
1003
# Line 997  class LinearPDE(object): Line 1016  class LinearPDE(object):
1016         sets a new solver package         sets a new solver package
1017
1018         @param package: sets a new solver method.         @param package: sets a new solver method.
1019         @type package: one of L{DEFAULT}, L{PASO} L{SCSL}, L{MKL}, L{UMFPACK}         @type package: one of L{DEFAULT}, L{PASO} L{SCSL}, L{MKL}, L{UMFPACK}, L{TRILINOS}
1020         """         """
1021         if package==None: package=self.DEFAULT         if package==None: package=self.DEFAULT
1022         if not package==self.getSolverPackage():         if not package==self.getSolverPackage():
# Line 1682  class LinearPDE(object): Line 1701  class LinearPDE(object):
1701                        raise ValueError,"coefficient B in lumped matrix may not be present."                        raise ValueError,"coefficient B in lumped matrix may not be present."
1702                   if not self.getCoefficientOfGeneralPDE("C").isEmpty():                   if not self.getCoefficientOfGeneralPDE("C").isEmpty():
1703                        raise ValueError,"coefficient C in lumped matrix may not be present."                        raise ValueError,"coefficient C in lumped matrix may not be present."
1704                     if not self.getCoefficientOfGeneralPDE("d_contact").isEmpty():
1705                          raise ValueError,"coefficient d_contact in lumped matrix may not be present."
1706                   if not self.getCoefficientOfGeneralPDE("A_reduced").isEmpty():                   if not self.getCoefficientOfGeneralPDE("A_reduced").isEmpty():
1707                        raise ValueError,"coefficient A_reduced in lumped matrix may not be present."                        raise ValueError,"coefficient A_reduced in lumped matrix may not be present."
1708                   if not self.getCoefficientOfGeneralPDE("B_reduced").isEmpty():                   if not self.getCoefficientOfGeneralPDE("B_reduced").isEmpty():
1709                        raise ValueError,"coefficient B_reduced in lumped matrix may not be present."                        raise ValueError,"coefficient B_reduced in lumped matrix may not be present."
1710                   if not self.getCoefficientOfGeneralPDE("C_reduced").isEmpty():                   if not self.getCoefficientOfGeneralPDE("C_reduced").isEmpty():
1711                        raise ValueError,"coefficient C_reduced in lumped matrix may not be present."                        raise ValueError,"coefficient C_reduced in lumped matrix may not be present."
1712                     if not self.getCoefficientOfGeneralPDE("d_contact_reduced").isEmpty():
1713                          raise ValueError,"coefficient d_contact_reduced in lumped matrix may not be present."
1714                   D=self.getCoefficientOfGeneralPDE("D")                   D=self.getCoefficientOfGeneralPDE("D")
1715                     d=self.getCoefficientOfGeneralPDE("d")
1716                     D_reduced=self.getCoefficientOfGeneralPDE("D_reduced")
1717                     d_reduced=self.getCoefficientOfGeneralPDE("d_reduced")
1718                   if not D.isEmpty():                   if not D.isEmpty():
1719                       if self.getNumSolutions()>1:                       if self.getNumSolutions()>1:
1720                          D_times_e=util.matrix_mult(D,numarray.ones((self.getNumSolutions(),)))                          D_times_e=util.matrix_mult(D,numarray.ones((self.getNumSolutions(),)))
# Line 1696  class LinearPDE(object): Line 1722  class LinearPDE(object):
1722                          D_times_e=D                          D_times_e=D
1723                   else:                   else:
1724                      D_times_e=escript.Data()                      D_times_e=escript.Data()
d=self.getCoefficientOfGeneralPDE("d")
1725                   if not d.isEmpty():                   if not d.isEmpty():
1726                       if self.getNumSolutions()>1:                       if self.getNumSolutions()>1:
1727                          d_times_e=util.matrix_mult(d,numarray.ones((self.getNumSolutions(),)))                          d_times_e=util.matrix_mult(d,numarray.ones((self.getNumSolutions(),)))
# Line 1704  class LinearPDE(object): Line 1729  class LinearPDE(object):
1729                          d_times_e=d                          d_times_e=d
1730                   else:                   else:
1731                      d_times_e=escript.Data()                      d_times_e=escript.Data()
1732                   d_contact=self.getCoefficientOfGeneralPDE("d_contact")
if not d_contact.isEmpty():
if self.getNumSolutions()>1:
d_contact_times_e=util.matrixmult(d_contact,numarray.ones((self.getNumSolutions(),)))
else:
d_contact_times_e=d_contact
else:
d_contact_times_e=escript.Data()

self.__operator=self.__getNewRightHandSide()
escript.Data(), \
D_times_e, \
d_times_e,\
d_contact_times_e)
D_reduced=self.getCoefficientOfGeneralPDE("D_reduced")
1733                   if not D_reduced.isEmpty():                   if not D_reduced.isEmpty():
1734                       if self.getNumSolutions()>1:                       if self.getNumSolutions()>1:
1735                          D_reduced_times_e=util.matrix_mult(D_reduced,numarray.ones((self.getNumSolutions(),)))                          D_reduced_times_e=util.matrix_mult(D_reduced,numarray.ones((self.getNumSolutions(),)))
# Line 1727  class LinearPDE(object): Line 1737  class LinearPDE(object):
1737                          D_reduced_times_e=D_reduced                          D_reduced_times_e=D_reduced
1738                   else:                   else:
1739                      D_reduced_times_e=escript.Data()                      D_reduced_times_e=escript.Data()
d_reduced=self.getCoefficientOfGeneralPDE("d_reduced")
1740                   if not d_reduced.isEmpty():                   if not d_reduced.isEmpty():
1741                       if self.getNumSolutions()>1:                       if self.getNumSolutions()>1:
1742                          d_reduced_times_e=util.matrix_mult(d_reduced,numarray.ones((self.getNumSolutions(),)))                          d_reduced_times_e=util.matrix_mult(d_reduced,numarray.ones((self.getNumSolutions(),)))
# Line 1735  class LinearPDE(object): Line 1744  class LinearPDE(object):
1744                          d_reduced_times_e=d_reduced                          d_reduced_times_e=d_reduced
1745                   else:                   else:
1746                      d_reduced_times_e=escript.Data()                      d_reduced_times_e=escript.Data()
1747                   d_contact_reduced=self.getCoefficientOfGeneralPDE("d_contact_reduced")
if not d_contact_reduced.isEmpty():
if self.getNumSolutions()>1:
d_contact_reduced_times_e=util.matrixmult(d_contact_reduced,numarray.ones((self.getNumSolutions(),)))
else:
d_contact_reduced_times_e=d_contact_reduced
else:
d_contact_reduced_times_e=escript.Data()

1748                   self.__operator=self.__getNewRightHandSide()                   self.__operator=self.__getNewRightHandSide()
1750                                                escript.Data(), \                      self.getDomain().addPDEToLumpedSystem(self.__operator, D_times_e, d_times_e)
1751                                                D_times_e, \                      self.getDomain().addPDEToLumpedSystem(self.__operator, D_reduced_times_e, d_reduced_times_e)
1752                                                d_times_e,\                   else:
1755                                                escript.Data(), \                                                   D_times_e, \
1756                                                D_reduced_times_e, \                                                   d_times_e,\
1757                                                d_reduced_times_e,\                                                   escript.Data())