/[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 150 by jgs, Thu Sep 15 03:44:45 2005 UTC revision 151 by jgs, Thu Sep 22 01:55:00 2005 UTC
# Line 2  Line 2 
2    
3  #  #
4  #      COPYRIGHT ACcESS 2004 -  All Rights Reserved  #      COPYRIGHT ACcESS 2004 -  All Rights Reserved
5  #  #
6  #   This software is the property of ACcESS.  No part of this code  #   This software is the property of ACcESS.  No part of this code
7  #   may be copied in any form or by any means without the expressed written  #   may be copied in any form or by any means without the expressed written
8  #   consent of ACcESS.  Copying, use or modification of this software  #   consent of ACcESS.  Copying, use or modification of this software
# Line 10  Line 10 
10  #   person has a software license agreement with ACcESS.  #   person has a software license agreement with ACcESS.
11  #  #
12  """  """
13  The module provides an interface to define and solve linear partial  The module provides an interface to define and solve linear partial
14  differential equations (PDEs) within L{escript}. L{linearPDEs} does not provide any  differential equations (PDEs) within L{escript}. L{linearPDEs} does not provide any
15  solver capabilities in itself but hands the PDE over to  solver capabilities in itself but hands the PDE over to
16  the PDE solver library defined through the L{Domain<escript.Domain>} of the PDE.  the PDE solver library defined through the L{Domain<escript.Domain>} of the PDE.
17  The general interface is provided through the L{LinearPDE} class. The  The general interface is provided through the L{LinearPDE} class. The
18  L{AdvectivePDE} which is derived from the L{LinearPDE} class  L{AdvectivePDE} which is derived from the L{LinearPDE} class
# Line 54  class UndefinedPDEError(ValueError): Line 54  class UndefinedPDEError(ValueError):
54     raised if a PDE is not fully defined yet.     raised if a PDE is not fully defined yet.
55     """     """
56    
57  class PDECoefficient:  class PDECoefficient(object):
58      """      """
59      A class for describing a PDE coefficient      A class for describing a PDE coefficient
60    
# Line 86  class PDECoefficient: Line 86  class PDECoefficient:
86      def __init__(self,where,pattern,altering):      def __init__(self,where,pattern,altering):
87         """         """
88         Initialise a PDE Coefficient type         Initialise a PDE Coefficient type
89          
90         @param where: describes where the coefficient lives         @param where: describes where the coefficient lives
91         @type where: one of L{INTERIOR}, L{BOUNDARY}, L{CONTACT}, L{SOLUTION}, L{REDUCED}         @type where: one of L{INTERIOR}, L{BOUNDARY}, L{CONTACT}, L{SOLUTION}, L{REDUCED}
92         @param pattern: describes the shape of the coefficient and how the shape is build for a given         @param pattern: describes the shape of the coefficient and how the shape is build for a given
93                spatial dimension and numbers of equation and solution in then PDE. For instance,                spatial dimension and numbers of equation and solution in then PDE. For instance,
94                (L{BY_EQUATION},L{BY_SOLUTION},L{BY_DIM}) descrbes a rank 3 coefficient which                (L{BY_EQUATION},L{BY_SOLUTION},L{BY_DIM}) descrbes a rank 3 coefficient which
95                is instanciated as shape (3,2,2) in case of a three equations and two solution components                is instanciated as shape (3,2,2) in case of a three equations and two solution components
# Line 101  class PDECoefficient: Line 101  class PDECoefficient:
101         @type altering: one of L{OPERATOR}, L{RIGHTHANDSIDE}, L{BOTH}         @type altering: one of L{OPERATOR}, L{RIGHTHANDSIDE}, L{BOTH}
102    
103         """         """
104           super(PDECoefficient, self).__init__()
105         self.what=where         self.what=where
106         self.pattern=pattern         self.pattern=pattern
107         self.altering=altering         self.altering=altering
# Line 125  class PDECoefficient: Line 126  class PDECoefficient:
126         @return:  L{FunctionSpace<escript.FunctionSpace>} of the coefficient         @return:  L{FunctionSpace<escript.FunctionSpace>} of the coefficient
127         @rtype:  L{FunctionSpace<escript.FunctionSpace>}         @rtype:  L{FunctionSpace<escript.FunctionSpace>}
128         """         """
129         if self.what==self.INTERIOR:         if self.what==self.INTERIOR:
130              return escript.Function(domain)              return escript.Function(domain)
131         elif self.what==self.BOUNDARY:         elif self.what==self.BOUNDARY:
132              return escript.FunctionOnBoundary(domain)              return escript.FunctionOnBoundary(domain)
133         elif self.what==self.CONTACT:         elif self.what==self.CONTACT:
134              return escript.FunctionOnContactZero(domain)              return escript.FunctionOnContactZero(domain)
135         elif self.what==self.SOLUTION:         elif self.what==self.SOLUTION:
136              if reducedEquationOrder and reducedSolutionOrder:              if reducedEquationOrder and reducedSolutionOrder:
137                  return escript.ReducedSolution(domain)                  return escript.ReducedSolution(domain)
138              else:              else:
139                  return escript.Solution(domain)                  return escript.Solution(domain)
140         elif self.what==self.REDUCED:         elif self.what==self.REDUCED:
141              if reducedEquationOrder and reducedSolutionOrder:              return escript.ReducedSolution(domain)
                 return escript.ReducedSolution(domain)  
             else:  
                 return escript.ReducedSolution(domain)  
142    
143      def getValue(self):      def getValue(self):
144         """         """
# Line 313  class PDECoefficient: Line 311  class PDECoefficient:
311                  s=s+(dim,)                  s=s+(dim,)
312         return s         return s
313    
314  class LinearPDE:  class LinearPDE(object):
315     """     """
316     This class is used to define a general linear, steady, second order PDE     This class is used to define a general linear, steady, second order PDE
317     for an unknown function M{u} on a given domain defined through a L{Domain<escript.Domain>} object.     for an unknown function M{u} on a given domain defined through a L{Domain<escript.Domain>} object.
318    
319     For a single PDE with a solution with a single component the linear PDE is defined in the following form:     For a single PDE with a solution with a single component the linear PDE is defined in the following form:
320      
321     M{-grad(A[j,l]*grad(u)[l]+B[j]u)[j]+C[l]*grad(u)[l]+D*u =-grad(X)[j,j]+Y}     M{-grad(A[j,l]*grad(u)[l]+B[j]u)[j]+C[l]*grad(u)[l]+D*u =-grad(X)[j,j]+Y}
322    
323     where M{grad(F)} denotes the spatial derivative of M{F}. Einstein's summation convention,     where M{grad(F)} denotes the spatial derivative of M{F}. Einstein's summation convention,
324     ie. summation over indexes appearing twice in a term of a sum is performed, is used.     ie. summation over indexes appearing twice in a term of a sum is performed, is used.
325     The coefficients M{A}, M{B}, M{C}, M{D}, M{X} and M{Y} have to be specified through L{Data<escript.Data>} objects in the     The coefficients M{A}, M{B}, M{C}, M{D}, M{X} and M{Y} have to be specified through L{Data<escript.Data>} objects in the
326     L{Function<escript.Function>} on the PDE or objects that can be converted into such L{Data<escript.Data>} objects.     L{Function<escript.Function>} on the PDE or objects that can be converted into such L{Data<escript.Data>} objects.
327     M{A} is a rank two, M{B}, M{C} and M{X} are rank one and M{D} and M{Y} are scalar.     M{A} is a rank two, M{B}, M{C} and M{X} are rank one and M{D} and M{Y} are scalar.
328    
329     The following natural boundary conditions are considered:     The following natural boundary conditions are considered:
330    
331     M{n[j]*(A[i,j]*grad(u)[l]+B[j]*u)+d*u=n[j]*X[j]+y}     M{n[j]*(A[i,j]*grad(u)[l]+B[j]*u)+d*u=n[j]*X[j]+y}
332    
333     where M{n} is the outer normal field calculated by L{getNormal<escript.FunctionSpace.getNormal>} of L{FunctionOnBoundary<escript.FunctionOnBoundary>}.     where M{n} is the outer normal field calculated by L{getNormal<escript.FunctionSpace.getNormal>} of L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
334     Notice that the coefficients M{A}, M{B} and M{X} are defined in the PDE. The coefficients M{d} and M{y} are       Notice that the coefficients M{A}, M{B} and M{X} are defined in the PDE. The coefficients M{d} and M{y} are
335     each a scalar in the L{FunctionOnBoundary<escript.FunctionOnBoundary>}.       each a scalar in the L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
336    
337    
338     Constraints for the solution prescribing the value of the solution at certain locations in the domain. They have the form     Constraints for the solution prescribing the value of the solution at certain locations in the domain. They have the form
# Line 343  class LinearPDE: Line 341  class LinearPDE:
341    
342     M{r} and M{q} are each scalar where M{q} is the characteristic function defining where the constraint is applied.     M{r} and M{q} are each scalar where M{q} is the characteristic function defining where the constraint is applied.
343     The constraints override any other condition set by the PDE or the boundary condition.     The constraints override any other condition set by the PDE or the boundary condition.
344      
345     The PDE is symmetrical if     The PDE is symmetrical if
346    
347     M{A[i,j]=A[j,i]}  and M{B[j]=C[j]}     M{A[i,j]=A[j,i]}  and M{B[j]=C[j]}
348    
349     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
350    
351     M{-grad(A[i,j,k,l]*grad(u[k])[l]+B[i,j,k]*u[k])[j]+C[i,k,l]*grad(u[k])[l]+D[i,k]*u[k] =-grad(X[i,j])[j]+Y[i] }     M{-grad(A[i,j,k,l]*grad(u[k])[l]+B[i,j,k]*u[k])[j]+C[i,k,l]*grad(u[k])[l]+D[i,k]*u[k] =-grad(X[i,j])[j]+Y[i] }
352    
353     M{A} is a ramk four, M{B} and M{C} are each a rank three, M{D} and M{X} are each a rank two and M{Y} is a rank one.     M{A} is a ramk four, M{B} and M{C} are each a rank three, M{D} and M{X} are each a rank two and M{Y} is a rank one.
354     The natural boundary conditions take the form:     The natural boundary conditions take the form:
355    
356     M{n[j]*(A[i,j,k,l]*grad(u[k])[l]+B[i,j,k]*u[k])+d[i,k]*u[k]=n[j]*X[i,j]+y[i]}     M{n[j]*(A[i,j,k,l]*grad(u[k])[l]+B[i,j,k]*u[k])+d[i,k]*u[k]=n[j]*X[i,j]+y[i]}
# Line 363  class LinearPDE: Line 361  class LinearPDE:
361    
362     M{u[i]=r[i]}  where  M{q[i]>0}     M{u[i]=r[i]}  where  M{q[i]>0}
363    
364     M{r} and M{q} are each rank one. Notice that at some locations not necessarily all components must have a constraint.     M{r} and M{q} are each rank one. Notice that at some locations not necessarily all components must have a constraint.
365    
366     The system of PDEs is symmetrical if     The system of PDEs is symmetrical if
367    
# Line 372  class LinearPDE: Line 370  class LinearPDE:
370          - M{D[i,k]=D[i,k]}          - M{D[i,k]=D[i,k]}
371          - M{d[i,k]=d[k,i]}          - M{d[i,k]=d[k,i]}
372    
373     L{LinearPDE} also supports solution discontinuities over a contact region in the domain. To specify the conditions across the     L{LinearPDE} also supports solution discontinuities over a contact region in the domain. To specify the conditions across the
374     discontinuity we are using the generalised flux M{J} which is in the case of a systems of PDEs and several components of the solution     discontinuity we are using the generalised flux M{J} which is in the case of a systems of PDEs and several components of the solution
375     defined as     defined as
376    
377     M{J[i,j]=A[i,j,k,l]*grad(u[k])[l]+B[i,j,k]*u[k]-X[i,j]}     M{J[i,j]=A[i,j,k,l]*grad(u[k])[l]+B[i,j,k]*u[k]-X[i,j]}
378    
# Line 387  class LinearPDE: Line 385  class LinearPDE:
385     the contact condition takes the form     the contact condition takes the form
386    
387     M{n[j]*J0[i,j]=n[j]*J1[i,j]=y_contact[i]- d_contact[i,k]*jump(u)[k]}     M{n[j]*J0[i,j]=n[j]*J1[i,j]=y_contact[i]- d_contact[i,k]*jump(u)[k]}
388      
389     where M{J0} and M{J1} are the fluxes on side 0 and side 1 of the discontinuity, respectively. M{jump(u)}, which is the difference     where M{J0} and M{J1} are the fluxes on side 0 and side 1 of the discontinuity, respectively. M{jump(u)}, which is the difference
390     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
391     L{jump<util.jump>}.     L{jump<util.jump>}.
392     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>}.
393     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
# Line 403  class LinearPDE: Line 401  class LinearPDE:
401     @cvar DIRECT: The direct solver based on LDU factorization     @cvar DIRECT: The direct solver based on LDU factorization
402     @cvar CHOLEVSKY: The direct solver based on LDLt factorization (can only be applied for symmetric PDEs)     @cvar CHOLEVSKY: The direct solver based on LDLt factorization (can only be applied for symmetric PDEs)
403     @cvar PCG: The preconditioned conjugate gradient method (can only be applied for symmetric PDEs)     @cvar PCG: The preconditioned conjugate gradient method (can only be applied for symmetric PDEs)
404     @cvar CR: The conjugate residual method     @cvar CR: The conjugate residual method
405     @cvar CGS: The conjugate gardient square method     @cvar CGS: The conjugate gardient square method
406     @cvar BICGSTAB: The stabilized BiConjugate Gradient method.     @cvar BICGSTAB: The stabilized BiConjugate Gradient method.
407     @cvar SSOR: The symmetric overrealaxtion method     @cvar SSOR: The symmetric overrealaxtion method
# Line 419  class LinearPDE: Line 417  class LinearPDE:
417     @cvar PASO: PASO solver package     @cvar PASO: PASO solver package
418     @cvar SCSL: SGI SCSL solver library     @cvar SCSL: SGI SCSL solver library
419     @cvar MKL: Intel's MKL solver library     @cvar MKL: Intel's MKL solver library
420     @cvar UMFPACK: the UMFPACK library     @cvar UMFPACK: the UMFPACK library
421     @cvar ITERATIVE: The default iterative solver     @cvar ITERATIVE: The default iterative solver
422    
423     """     """
# Line 466  class LinearPDE: Line 464  class LinearPDE:
464       @param debug: if True debug informations are printed.       @param debug: if True debug informations are printed.
465    
466       """       """
467         super(LinearPDE, self).__init__()
468       #       #
469       #   the coefficients of the general PDE:       #   the coefficients of the general PDE:
470       #       #
# Line 602  class LinearPDE: Line 601  class LinearPDE:
601       @rtype: L{bool}       @rtype: L{bool}
602       """       """
603       return self.__reduce_solution_order       return self.__reduce_solution_order
604    
605     def getFunctionSpaceForEquation(self):     def getFunctionSpaceForEquation(self):
606       """       """
607       returns the L{FunctionSpace<escript.FunctionSpace>} used to discretize the equation       returns the L{FunctionSpace<escript.FunctionSpace>} used to discretize the equation
# Line 802  class LinearPDE: Line 801  class LinearPDE:
801    
802       M{J[i,j]=A[i,j,k,l]*grad(u[k])[l]+B[i,j,k]u[k]-X[i,j]}       M{J[i,j]=A[i,j,k,l]*grad(u[k])[l]+B[i,j,k]u[k]-X[i,j]}
803    
804       or       or
805    
806       M{J[j]=A[i,j]*grad(u)[l]+B[j]u-X[j]}       M{J[j]=A[i,j]*grad(u)[l]+B[j]u-X[j]}
807    
# Line 865  class LinearPDE: Line 864  class LinearPDE:
864         """         """
865         returns the solver method         returns the solver method
866    
867         @return: the solver method currently be used.         @return: the solver method currently be used.
868         @rtype: C{int}         @rtype: C{int}
869         """         """
870         return self.__solver_method         return self.__solver_method
# Line 887  class LinearPDE: Line 886  class LinearPDE:
886         """         """
887         returns the package of the solver         returns the package of the solver
888    
889         @return: the solver package currently being used.         @return: the solver package currently being used.
890         @rtype: C{int}         @rtype: C{int}
891         """         """
892         return self.__solver_package         return self.__solver_package
# Line 1614  class Poisson(LinearPDE): Line 1613  class Poisson(LinearPDE):
1613       @param debug: if True debug informations are printed.       @param debug: if True debug informations are printed.
1614    
1615       """       """
1616       LinearPDE.__init__(self,domain,1,1,debug)       super(Poisson, self).__init__(domain,1,1,debug)
1617       self.COEFFICIENTS={"f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),       self.COEFFICIENTS={"f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1618                            "q": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}                            "q": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}
1619       self.setSymmetryOn()       self.setSymmetryOn()
# Line 1631  class Poisson(LinearPDE): Line 1630  class Poisson(LinearPDE):
1630                 depending of reduced order is used for the representation of the equation.                 depending of reduced order is used for the representation of the equation.
1631       @raise IllegalCoefficient: if an unknown coefficient keyword is used.       @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1632       """       """
1633       LinearPDE.setValue(self,**coefficients)       super(Poisson, self).setValue(**coefficients)
1634    
1635     def getCoefficientOfGeneralPDE(self,name):     def getCoefficientOfGeneralPDE(self,name):
1636       """       """
# Line 1696  class Helmholtz(LinearPDE): Line 1695  class Helmholtz(LinearPDE):
1695       @param debug: if True debug informations are printed.       @param debug: if True debug informations are printed.
1696    
1697       """       """
1698       LinearPDE.__init__(self,domain,1,1,debug)       super(Helmholtz, self).__init__(domain,1,1,debug)
1699       self.COEFFICIENTS={"omega": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),       self.COEFFICIENTS={"omega": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
1700                          "k": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),                          "k": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
1701                          "f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),                          "f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
# Line 1729  class Helmholtz(LinearPDE): Line 1728  class Helmholtz(LinearPDE):
1728                 depending of reduced order is used for the representation of the equation.                 depending of reduced order is used for the representation of the equation.
1729       @raise IllegalCoefficient: if an unknown coefficient keyword is used.       @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1730       """       """
1731       LinearPDE.setValue(self,**coefficients)       super(Helmholtz, self).setValue(**coefficients)
1732    
1733     def getCoefficientOfGeneralPDE(self,name):     def getCoefficientOfGeneralPDE(self,name):
1734       """       """
# Line 1787  class LameEquation(LinearPDE): Line 1786  class LameEquation(LinearPDE):
1786     """     """
1787    
1788     def __init__(self,domain,debug=False):     def __init__(self,domain,debug=False):
1789         LinearPDE.__init__(self,domain,domain.getDim(),domain.getDim(),debug)        super(LameEquation, self).__init__(domain,\
1790         self.COEFFICIENTS={ "lame_lambda"  : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),                                           domain.getDim(),domain.getDim(),debug)
1791          self.COEFFICIENTS={ "lame_lambda"  : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),
1792                            "lame_mu"      : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),                            "lame_mu"      : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),
1793                            "F"            : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),                            "F"            : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1794                            "sigma"        : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM),PDECoefficient.RIGHTHANDSIDE),                            "sigma"        : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM),PDECoefficient.RIGHTHANDSIDE),
1795                            "f"            : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),                            "f"            : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1796                            "r"            : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH),                            "r"            : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH),
1797                            "q"            : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}                            "q"            : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}
1798         self.setSymmetryOn()        self.setSymmetryOn()
1799    
1800     def setValue(self,**coefficients):     def setValue(self,**coefficients):
1801       """       """
# Line 1820  class LameEquation(LinearPDE): Line 1820  class LameEquation(LinearPDE):
1820                 depending of reduced order is used for the representation of the equation.                 depending of reduced order is used for the representation of the equation.
1821       @raise IllegalCoefficient: if an unknown coefficient keyword is used.       @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1822       """       """
1823       LinearPDE.setValue(self,**coefficients)       super(LameEquation, self).setValue(**coefficients)
1824    
1825     def getCoefficientOfGeneralPDE(self,name):     def getCoefficientOfGeneralPDE(self,name):
1826       """       """
# Line 1876  class AdvectivePDE(LinearPDE): Line 1876  class AdvectivePDE(LinearPDE):
1876    
1877     M{Z[j]=C[j]-B[j]}     M{Z[j]=C[j]-B[j]}
1878    
1879     or     or
1880    
1881     M{Z[i,k,l]=C[i,k,l]-B[i,l,k]}     M{Z[i,k,l]=C[i,k,l]-B[i,l,k]}
1882    
1883     To measure the dominance of the advective terms over the diffusive term M{A} the     To measure the dominance of the advective terms over the diffusive term M{A} the
1884     X{Pelclet number} M{P} is used. It is defined as     X{Pelclet number} M{P} is used. It is defined as
1885    
1886     M{P=h|Z|/(2|A|)}     M{P=h|Z|/(2|A|)}
1887    
1888     where M{|.|} denotes the L{length<util.length>} of the arument and M{h} is the local cell size     where M{|.|} denotes the L{length<util.length>} of the arument and M{h} is the local cell size
1889     from L{getSize<escript.Domain.getSize>}. Where M{|A|==0} M{P} is M{S{infinity}}.     from L{getSize<escript.Domain.getSize>}. Where M{|A|==0} M{P} is M{S{infinity}}.
1890    
1891     From the X{Pelclet number} the stabilization parameters M{S{Xi}} and M{S{Xi}} are calculated:     From the X{Pelclet number} the stabilization parameters M{S{Xi}} and M{S{Xi}} are calculated:
# Line 1913  class AdvectivePDE(LinearPDE): Line 1913  class AdvectivePDE(LinearPDE):
1913     """     """
1914     def __init__(self,domain,numEquations=None,numSolutions=None,xi=None,debug=False):     def __init__(self,domain,numEquations=None,numSolutions=None,xi=None,debug=False):
1915        """        """
1916        creates a linear, steady, second order PDE on a L{Domain<escript.Domain>}        creates a linear, steady, second order PDE on a L{Domain<escript.Domain>}
1917    
1918        @param domain: domain of the PDE        @param domain: domain of the PDE
1919        @type domain: L{Domain<escript.Domain>}        @type domain: L{Domain<escript.Domain>}
# Line 1921  class AdvectivePDE(LinearPDE): Line 1921  class AdvectivePDE(LinearPDE):
1921                             is exracted from the PDE coefficients.                             is exracted from the PDE coefficients.
1922        @param numSolutions: number of solution components. If  numSolutions==None the number of solution components        @param numSolutions: number of solution components. If  numSolutions==None the number of solution components
1923                             is exracted from the PDE coefficients.                             is exracted from the PDE coefficients.
1924        @param xi: defines a function which returns for any given Preclet number as L{Scalar<escript.Scalar>} object the        @param xi: defines a function which returns for any given Preclet number as L{Scalar<escript.Scalar>} object the
1925                   M{S{xi}}-value used to define the stabilization parameters. If equal to None, L{ELMAN_RAMAGE} is used.                   M{S{xi}}-value used to define the stabilization parameters. If equal to None, L{ELMAN_RAMAGE} is used.
1926        @type xi: callable object which returns a L{Scalar<escript.Scalar>} object.        @type xi: callable object which returns a L{Scalar<escript.Scalar>} object.
1927        @param debug: if True debug informations are printed.        @param debug: if True debug informations are printed.
1928        """        """
1929          super(AdvectivePDE, self).__init__(domain,\
1930        LinearPDE.__init__(self,domain,numEquations,numSolutions,debug)                                           numEquations,numSolutions,debug)
1931        if xi==None:        if xi==None:
1932           self.__xi=AdvectivePDE.ELMAN_RAMAGE           self.__xi=AdvectivePDE.ELMAN_RAMAGE
1933        else:        else:
1934           self.__xi=xi           self.__xi=xi
1935        self.__Xi=escript.Data()        self.__Xi=escript.Data()
1936    
1937     def setValue(self,**coefficients):     def setValue(**coefficients):
1938        """        """
1939        sets new values to coefficients        sets new values to coefficients
1940    
# Line 1971  class AdvectivePDE(LinearPDE): Line 1971  class AdvectivePDE(LinearPDE):
1971    
1972        """        """
1973        if "A" in coefficients.keys()   or "B" in coefficients.keys() or "C" in coefficients.keys(): self.__Xi=escript.Data()        if "A" in coefficients.keys()   or "B" in coefficients.keys() or "C" in coefficients.keys(): self.__Xi=escript.Data()
1974        LinearPDE.setValue(self,**coefficients)        super(AdvectivePDE, self).setValue(**coefficients)
1975    
1976     def ELMAN_RAMAGE(self,P):     def ELMAN_RAMAGE(self,P):
1977       """       """
1978       Predefined function to set a values for M{S{xi}} from a Preclet number M{P}.       Predefined function to set a values for M{S{xi}} from a Preclet number M{P}.
1979       This function uses the method suggested by H.C. Elman and A. Ramage, I{SIAM J. Numer. Anal.}, B{40} (2002)       This function uses the method suggested by H.C. Elman and A. Ramage, I{SIAM J. Numer. Anal.}, B{40} (2002)
1980            - M{S{xi}(P)=0} for M{P<1}            - M{S{xi}(P)=0} for M{P<1}
1981            - M{S{xi}(P)=(1-1/P)/2} otherwise            - M{S{xi}(P)=(1-1/P)/2} otherwise
1982    
1983       @param P: Preclet number       @param P: Preclet number
1984       @type P: L{Scalar<escript.Scalar>}       @type P: L{Scalar<escript.Scalar>}
1985       @return: up-wind weightimg factor       @return: up-wind weightimg factor
1986       @rtype: L{Scalar<escript.Scalar>}       @rtype: L{Scalar<escript.Scalar>}
1987       """       """
1988       return (P-1.).wherePositive()*0.5*(1.-1./(P+1.e-15))       return (P-1.).wherePositive()*0.5*(1.-1./(P+1.e-15))
1989    
1990     def SIMPLIFIED_BROOK_HUGHES(self,P):     def SIMPLIFIED_BROOK_HUGHES(self,P):
1991       """       """
1992       Predefined function to set a values for M{S{xi}} from a Preclet number M{P}.       Predefined function to set a values for M{S{xi}} from a Preclet number M{P}.
1993       The original methods is       The original methods is
1994        
1995       M{S{xi}(P)=coth(P)-1/P}       M{S{xi}(P)=coth(P)-1/P}
1996    
1997       As the evaluation of M{coth} is expensive we are using the approximation:       As the evaluation of M{coth} is expensive we are using the approximation:
1998        
1999           - M{S{xi}(P)=P/3} where M{P<3}           - M{S{xi}(P)=P/3} where M{P<3}
2000           - M{S{xi}(P)=1/2} otherwise           - M{S{xi}(P)=1/2} otherwise
2001    
2002       @param P: Preclet number       @param P: Preclet number
2003       @type P: L{Scalar<escript.Scalar>}       @type P: L{Scalar<escript.Scalar>}
2004       @return: up-wind weightimg factor       @return: up-wind weightimg factor
2005       @rtype: L{Scalar<escript.Scalar>}       @rtype: L{Scalar<escript.Scalar>}
2006       """       """
2007       c=(P-3.).whereNegative()       c=(P-3.).whereNegative()
2008       return P/6.*c+1./2.*(1.-c)       return P/6.*c+1./2.*(1.-c)
2009    
2010     def HALF(self,P):     def HALF(self,P):
2011       """       """
2012       Predefined function to set value M{1/2} for M{S{xi}}       Predefined function to set value M{1/2} for M{S{xi}}
2013    
2014       @param P: Preclet number       @param P: Preclet number
2015       @type P: L{Scalar<escript.Scalar>}       @type P: L{Scalar<escript.Scalar>}
2016       @return: up-wind weightimg factor       @return: up-wind weightimg factor
2017       @rtype: L{Scalar<escript.Scalar>}       @rtype: L{Scalar<escript.Scalar>}
2018       """       """
2019       return escript.Scalar(0.5,P.getFunctionSpace())       return escript.Scalar(0.5,P.getFunctionSpace())
2020    
2021     def __calculateXi(self,peclet_factor,Z,h):     def __calculateXi(self,peclet_factor,Z,h):
# Line 2219  class AdvectivePDE(LinearPDE): Line 2219  class AdvectivePDE(LinearPDE):
2219    
2220    
2221  # $Log$  # $Log$
2222    # Revision 1.14  2005/09/22 01:54:57  jgs
2223    # Merge of development branch dev-02 back to main trunk on 2005-09-22
2224    #
2225  # Revision 1.13  2005/09/15 03:44:19  jgs  # Revision 1.13  2005/09/15 03:44:19  jgs
2226  # Merge of development branch dev-02 back to main trunk on 2005-09-15  # Merge of development branch dev-02 back to main trunk on 2005-09-15
2227  #  #
# Line 2231  class AdvectivePDE(LinearPDE): Line 2234  class AdvectivePDE(LinearPDE):
2234  # Revision 1.10  2005/08/12 01:45:36  jgs  # Revision 1.10  2005/08/12 01:45:36  jgs
2235  # erge of development branch dev-02 back to main trunk on 2005-08-12  # erge of development branch dev-02 back to main trunk on 2005-08-12
2236  #  #
2237    # Revision 1.9.2.17  2005/09/21 07:03:33  matt
2238    # PDECoefficient and LinearPDE are now new style classes (introduced in Python
2239    # 2.2). Classes Poisson, Helmholtz, LameEquation and AdvectivePDE have been
2240    # modified to instead use portable/cooperative "super" calls to extend base
2241    # class methods.
2242    #
2243    # Revision 1.9.2.16  2005/09/16 01:54:37  matt
2244    # Removed redundant if-loop.
2245    #
2246  # Revision 1.9.2.15  2005/09/14 08:09:18  matt  # Revision 1.9.2.15  2005/09/14 08:09:18  matt
2247  # Added "REDUCED" solution PDECoefficient descriptors for LinearPDEs.  # Added "REDUCED" solution PDECoefficient descriptors for LinearPDEs.
2248  #  #

Legend:
Removed from v.150  
changed lines
  Added in v.151

  ViewVC Help
Powered by ViewVC 1.1.26