/[escript]/temp_trunk_copy/escript/py_src/linearPDEs.py
ViewVC logotype

Diff of /temp_trunk_copy/escript/py_src/linearPDEs.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/esys2/escript/py_src/linearPDEs.py revision 150 by jgs, Thu Sep 15 03:44:45 2005 UTC trunk/escript/py_src/linearPDEs.py revision 430 by gross, Wed Jan 11 06:40:50 2006 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
19  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},
20  L{Helmholtz}, L{LameEquation}  L{Helmholtz}, L{LameEquation}, L{AdvectionDiffusion}
21  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
22  to define of solve these sepecial PDEs.  to define of solve these sepecial PDEs.
23    
# 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       @cvar AMG: algebraic multi grid
423       @cvar RILU: recursive ILU
424    
425     """     """
426     DEFAULT= 0     DEFAULT= 0
# Line 445  class LinearPDE: Line 445  class LinearPDE:
445     UMFPACK= 16     UMFPACK= 16
446     ITERATIVE= 20     ITERATIVE= 20
447     PASO= 21     PASO= 21
448       AMG= 22
449       RILU = 23
450    
451     __TOL=1.e-13     __TOL=1.e-13
452     __PACKAGE_KEY="package"     __PACKAGE_KEY="package"
453     __METHOD_KEY="method"     __METHOD_KEY="method"
454     __SYMMETRY_KEY="symmetric"     __SYMMETRY_KEY="symmetric"
455     __TOLERANCE_KEY="tolerance"     __TOLERANCE_KEY="tolerance"
456       __PRECONDITIONER_KEY="preconditioner"
457    
458    
459     def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):     def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):
# Line 466  class LinearPDE: Line 469  class LinearPDE:
469       @param debug: if True debug informations are printed.       @param debug: if True debug informations are printed.
470    
471       """       """
472         super(LinearPDE, self).__init__()
473       #       #
474       #   the coefficients of the general PDE:       #   the coefficients of the general PDE:
475       #       #
# Line 499  class LinearPDE: Line 503  class LinearPDE:
503       self.__tolerance=1.e-8       self.__tolerance=1.e-8
504       self.__solver_method=self.DEFAULT       self.__solver_method=self.DEFAULT
505       self.__solver_package=self.DEFAULT       self.__solver_package=self.DEFAULT
506         self.__preconditioner=self.DEFAULT
507       self.__matrix_type=self.__domain.getSystemMatrixTypeId(self.DEFAULT,self.DEFAULT,False)       self.__matrix_type=self.__domain.getSystemMatrixTypeId(self.DEFAULT,self.DEFAULT,False)
508       self.__sym=False       self.__sym=False
509    
# Line 602  class LinearPDE: Line 607  class LinearPDE:
607       @rtype: L{bool}       @rtype: L{bool}
608       """       """
609       return self.__reduce_solution_order       return self.__reduce_solution_order
610    
611     def getFunctionSpaceForEquation(self):     def getFunctionSpaceForEquation(self):
612       """       """
613       returns the L{FunctionSpace<escript.FunctionSpace>} used to discretize the equation       returns the L{FunctionSpace<escript.FunctionSpace>} used to discretize the equation
# Line 773  class LinearPDE: Line 778  class LinearPDE:
778         @type verbose: C{bool}         @type verbose: C{bool}
779         @keyword reordering: reordering scheme to be used during elimination. Allowed values are         @keyword reordering: reordering scheme to be used during elimination. Allowed values are
780                              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}  
781         @keyword iter_max: maximum number of iteration steps allowed.         @keyword iter_max: maximum number of iteration steps allowed.
782         @keyword drop_tolerance: threshold for drupping in L{ILUT}         @keyword drop_tolerance: threshold for drupping in L{ILUT}
783         @keyword drop_storage: maximum of allowed memory in L{ILUT}         @keyword drop_storage: maximum of allowed memory in L{ILUT}
# Line 787  class LinearPDE: Line 790  class LinearPDE:
790               self.__solution=self.copyConstraint(f*mat)               self.__solution=self.copyConstraint(f*mat)
791            else:            else:
792               options[self.__TOLERANCE_KEY]=self.getTolerance()               options[self.__TOLERANCE_KEY]=self.getTolerance()
793               options[self.__METHOD_KEY]=self.getSolverMethod()               options[self.__METHOD_KEY]=self.getSolverMethod()[0]
794                 options[self.__PRECONDITIONER_KEY]=self.getSolverMethod()[1]
795               options[self.__PACKAGE_KEY]=self.getSolverPackage()               options[self.__PACKAGE_KEY]=self.getSolverPackage()
796               options[self.__SYMMETRY_KEY]=self.isSymmetric()               options[self.__SYMMETRY_KEY]=self.isSymmetric()
797               self.trace("PDE is resolved.")               self.trace("PDE is resolved.")
# Line 802  class LinearPDE: Line 806  class LinearPDE:
806    
807       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]}
808    
809       or       or
810    
811       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]}
812    
# Line 816  class LinearPDE: Line 820  class LinearPDE:
820     # =============================================================================     # =============================================================================
821     #   solver settings:     #   solver settings:
822     # =============================================================================     # =============================================================================
823     def setSolverMethod(self,solver=None):     def setSolverMethod(self,solver=None,preconditioner=None):
824         """         """
825         sets a new solver         sets a new solver
826    
827         @param solver: sets a new solver method.         @param solver: sets a new solver method.
828         @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}.
829           @param preconditioner: sets a new solver method.
830           @type solver: one of L{DEFAULT}, L{JACOBI} L{ILU0}, L{ILUT},L{SSOR}
831         """         """
832         if solver==None: solve=self.DEFAULT         if solver==None: solve=self.DEFAULT
833         if not solver==self.getSolverMethod():         if preconditioner==None: preconditioner=self.DEFAULT
834           if not (solver,preconditioner)==self.getSolverMethod():
835             self.__solver_method=solver             self.__solver_method=solver
836               self.__preconditioner=preconditioner
837             self.__checkMatrixType()             self.__checkMatrixType()
838             self.trace("New solver is %s"%self.getSolverMethodName())             self.trace("New solver is %s"%self.getSolverMethodName())
839    
# Line 839  class LinearPDE: Line 847  class LinearPDE:
847    
848         m=self.getSolverMethod()         m=self.getSolverMethod()
849         p=self.getSolverPackage()         p=self.getSolverPackage()
850         if m==self.DEFAULT: method="DEFAULT"         method=""
851         elif m==self.DIRECT: method= "DIRECT"         if m[0]==self.DEFAULT: method="DEFAULT"
852         elif m==self.ITERATIVE: method= "ITERATIVE"         elif m[0]==self.DIRECT: method= "DIRECT"
853         elif m==self.CHOLEVSKY: method= "CHOLEVSKY"         elif m[0]==self.ITERATIVE: method= "ITERATIVE"
854         elif m==self.PCG: method= "PCG"         elif m[0]==self.CHOLEVSKY: method= "CHOLEVSKY"
855         elif m==self.CR: method= "CR"         elif m[0]==self.PCG: method= "PCG"
856         elif m==self.CGS: method= "CGS"         elif m[0]==self.CR: method= "CR"
857         elif m==self.BICGSTAB: method= "BICGSTAB"         elif m[0]==self.CGS: method= "CGS"
858         elif m==self.SSOR: method= "SSOR"         elif m[0]==self.BICGSTAB: method= "BICGSTAB"
859         elif m==self.GMRES: method= "GMRES"         elif m[0]==self.SSOR: method= "SSOR"
860         elif m==self.PRES20: method= "PRES20"         elif m[0]==self.GMRES: method= "GMRES"
861         elif m==self.LUMPING: method= "LUMPING"         elif m[0]==self.PRES20: method= "PRES20"
862         else : method="unknown"         elif m[0]==self.LUMPING: method= "LUMPING"
863           if m[1]==self.DEFAULT: method+="+DEFAULT"
864           elif m[1]==self.JACOBI: method+= "+JACOBI"
865           elif m[1]==self.ILU0: method+= "+ILU0"
866           elif m[1]==self.ILUT: method+= "+ILUT"
867           elif m[1]==self.SSOR: method+= "+SSOR"
868         if p==self.DEFAULT: package="DEFAULT"         if p==self.DEFAULT: package="DEFAULT"
869         elif p==self.PASO: package= "PASO"         elif p==self.PASO: package= "PASO"
870         elif p==self.MKL: package= "MKL"         elif p==self.MKL: package= "MKL"
# Line 865  class LinearPDE: Line 878  class LinearPDE:
878         """         """
879         returns the solver method         returns the solver method
880    
881         @return: the solver method currently be used.         @return: the solver method currently be used.
882         @rtype: C{int}         @rtype: C{int}
883         """         """
884         return self.__solver_method         return self.__solver_method,self.__preconditioner
885    
886     def setSolverPackage(self,package=None):     def setSolverPackage(self,package=None):
887         """         """
# Line 887  class LinearPDE: Line 900  class LinearPDE:
900         """         """
901         returns the package of the solver         returns the package of the solver
902    
903         @return: the solver package currently being used.         @return: the solver package currently being used.
904         @rtype: C{int}         @rtype: C{int}
905         """         """
906         return self.__solver_package         return self.__solver_package
# Line 899  class LinearPDE: Line 912  class LinearPDE:
912        @return: True is lumping is currently used a solver method.        @return: True is lumping is currently used a solver method.
913        @rtype: C{bool}        @rtype: C{bool}
914        """        """
915        return self.getSolverMethod()==self.LUMPING        return self.getSolverMethod()[0]==self.LUMPING
916    
917     def setTolerance(self,tol=1.e-8):     def setTolerance(self,tol=1.e-8):
918         """         """
# Line 1092  class LinearPDE: Line 1105  class LinearPDE:
1105       """       """
1106       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.
1107       """       """
1108       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())
1109       if not new_matrix_type==self.__matrix_type:       if not new_matrix_type==self.__matrix_type:
1110           self.trace("Matrix type is now %d."%new_matrix_type)           self.trace("Matrix type is now %d."%new_matrix_type)
1111           self.__matrix_type=new_matrix_type           self.__matrix_type=new_matrix_type
# Line 1605  class Poisson(LinearPDE): Line 1618  class Poisson(LinearPDE):
1618    
1619     """     """
1620    
1621     def __init__(self,domain,f=escript.Data(),q=escript.Data(),debug=False):     def __init__(self,domain,debug=False):
1622       """       """
1623       initializes a new Poisson equation       initializes a new Poisson equation
1624    
# Line 1614  class Poisson(LinearPDE): Line 1627  class Poisson(LinearPDE):
1627       @param debug: if True debug informations are printed.       @param debug: if True debug informations are printed.
1628    
1629       """       """
1630       LinearPDE.__init__(self,domain,1,1,debug)       super(Poisson, self).__init__(domain,1,1,debug)
1631       self.COEFFICIENTS={"f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),       self.COEFFICIENTS={"f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1632                            "q": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}                            "q": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}
1633       self.setSymmetryOn()       self.setSymmetryOn()
# Line 1631  class Poisson(LinearPDE): Line 1644  class Poisson(LinearPDE):
1644                 depending of reduced order is used for the representation of the equation.                 depending of reduced order is used for the representation of the equation.
1645       @raise IllegalCoefficient: if an unknown coefficient keyword is used.       @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1646       """       """
1647       LinearPDE.setValue(self,**coefficients)       super(Poisson, self).setValue(**coefficients)
1648    
1649     def getCoefficientOfGeneralPDE(self,name):     def getCoefficientOfGeneralPDE(self,name):
1650       """       """
# Line 1645  class Poisson(LinearPDE): Line 1658  class Poisson(LinearPDE):
1658       @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.       @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
1659       """       """
1660       if name == "A" :       if name == "A" :
1661           return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))           return escript.Data(util.kronecker(self.getDim()),escript.Function(self.getDomain()))
1662       elif name == "B" :       elif name == "B" :
1663           return escript.Data()           return escript.Data()
1664       elif name == "C" :       elif name == "C" :
# Line 1696  class Helmholtz(LinearPDE): Line 1709  class Helmholtz(LinearPDE):
1709       @param debug: if True debug informations are printed.       @param debug: if True debug informations are printed.
1710    
1711       """       """
1712       LinearPDE.__init__(self,domain,1,1,debug)       super(Helmholtz, self).__init__(domain,1,1,debug)
1713       self.COEFFICIENTS={"omega": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),       self.COEFFICIENTS={"omega": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
1714                          "k": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),                          "k": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
1715                          "f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),                          "f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
# Line 1729  class Helmholtz(LinearPDE): Line 1742  class Helmholtz(LinearPDE):
1742                 depending of reduced order is used for the representation of the equation.                 depending of reduced order is used for the representation of the equation.
1743       @raise IllegalCoefficient: if an unknown coefficient keyword is used.       @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1744       """       """
1745       LinearPDE.setValue(self,**coefficients)       super(Helmholtz, self).setValue(**coefficients)
1746    
1747     def getCoefficientOfGeneralPDE(self,name):     def getCoefficientOfGeneralPDE(self,name):
1748       """       """
# Line 1787  class LameEquation(LinearPDE): Line 1800  class LameEquation(LinearPDE):
1800     """     """
1801    
1802     def __init__(self,domain,debug=False):     def __init__(self,domain,debug=False):
1803         LinearPDE.__init__(self,domain,domain.getDim(),domain.getDim(),debug)        super(LameEquation, self).__init__(domain,\
1804         self.COEFFICIENTS={ "lame_lambda"  : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),                                           domain.getDim(),domain.getDim(),debug)
1805          self.COEFFICIENTS={ "lame_lambda"  : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),
1806                            "lame_mu"      : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),                            "lame_mu"      : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),
1807                            "F"            : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),                            "F"            : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1808                            "sigma"        : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM),PDECoefficient.RIGHTHANDSIDE),                            "sigma"        : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM),PDECoefficient.RIGHTHANDSIDE),
1809                            "f"            : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),                            "f"            : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1810                            "r"            : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH),                            "r"            : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH),
1811                            "q"            : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}                            "q"            : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}
1812         self.setSymmetryOn()        self.setSymmetryOn()
1813    
1814     def setValue(self,**coefficients):     def setValue(self,**coefficients):
1815       """       """
# Line 1820  class LameEquation(LinearPDE): Line 1834  class LameEquation(LinearPDE):
1834                 depending of reduced order is used for the representation of the equation.                 depending of reduced order is used for the representation of the equation.
1835       @raise IllegalCoefficient: if an unknown coefficient keyword is used.       @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1836       """       """
1837       LinearPDE.setValue(self,**coefficients)       super(LameEquation, self).setValue(**coefficients)
1838    
1839     def getCoefficientOfGeneralPDE(self,name):     def getCoefficientOfGeneralPDE(self,name):
1840       """       """
# Line 1876  class AdvectivePDE(LinearPDE): Line 1890  class AdvectivePDE(LinearPDE):
1890    
1891     M{Z[j]=C[j]-B[j]}     M{Z[j]=C[j]-B[j]}
1892    
1893     or     or
1894    
1895     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]}
1896    
1897     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
1898     X{Pelclet number} M{P} is used. It is defined as     X{Pelclet number} M{P} is used. It is defined as
1899    
1900     M{P=h|Z|/(2|A|)}     M{P=h|Z|/(2|A|)}
1901    
1902     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
1903     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}}.
1904    
1905     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 1927  class AdvectivePDE(LinearPDE):
1927     """     """
1928     def __init__(self,domain,numEquations=None,numSolutions=None,xi=None,debug=False):     def __init__(self,domain,numEquations=None,numSolutions=None,xi=None,debug=False):
1929        """        """
1930        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>}
1931    
1932        @param domain: domain of the PDE        @param domain: domain of the PDE
1933        @type domain: L{Domain<escript.Domain>}        @type domain: L{Domain<escript.Domain>}
# Line 1921  class AdvectivePDE(LinearPDE): Line 1935  class AdvectivePDE(LinearPDE):
1935                             is exracted from the PDE coefficients.                             is exracted from the PDE coefficients.
1936        @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
1937                             is exracted from the PDE coefficients.                             is exracted from the PDE coefficients.
1938        @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
1939                   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.
1940        @type xi: callable object which returns a L{Scalar<escript.Scalar>} object.        @type xi: callable object which returns a L{Scalar<escript.Scalar>} object.
1941        @param debug: if True debug informations are printed.        @param debug: if True debug informations are printed.
1942        """        """
1943          super(AdvectivePDE, self).__init__(domain,\
1944        LinearPDE.__init__(self,domain,numEquations,numSolutions,debug)                                           numEquations,numSolutions,debug)
1945        if xi==None:        if xi==None:
1946           self.__xi=AdvectivePDE.ELMAN_RAMAGE           self.__xi=AdvectivePDE.ELMAN_RAMAGE
1947        else:        else:
1948           self.__xi=xi           self.__xi=xi
1949        self.__Xi=escript.Data()        self.__Xi=escript.Data()
1950    
1951     def setValue(self,**coefficients):     def setValue(**coefficients):
1952        """        """
1953        sets new values to coefficients        sets new values to coefficients
1954    
# Line 1971  class AdvectivePDE(LinearPDE): Line 1985  class AdvectivePDE(LinearPDE):
1985    
1986        """        """
1987        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()
1988        LinearPDE.setValue(self,**coefficients)        super(AdvectivePDE, self).setValue(**coefficients)
1989    
1990     def ELMAN_RAMAGE(self,P):     def ELMAN_RAMAGE(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       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)
1994            - M{S{xi}(P)=0} for M{P<1}            - M{S{xi}(P)=0} for M{P<1}
1995            - M{S{xi}(P)=(1-1/P)/2} otherwise            - M{S{xi}(P)=(1-1/P)/2} otherwise
1996    
1997       @param P: Preclet number       @param P: Preclet number
1998       @type P: L{Scalar<escript.Scalar>}       @type P: L{Scalar<escript.Scalar>}
1999       @return: up-wind weightimg factor       @return: up-wind weightimg factor
2000       @rtype: L{Scalar<escript.Scalar>}       @rtype: L{Scalar<escript.Scalar>}
2001       """       """
2002       return (P-1.).wherePositive()*0.5*(1.-1./(P+1.e-15))       return util.wherePositive(P-1.)*0.5*(1.-1./(P+1.e-15))
2003    
2004     def SIMPLIFIED_BROOK_HUGHES(self,P):     def SIMPLIFIED_BROOK_HUGHES(self,P):
2005       """       """
2006       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}.
2007       The original methods is       The original methods is
2008        
2009       M{S{xi}(P)=coth(P)-1/P}       M{S{xi}(P)=coth(P)-1/P}
2010    
2011       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:
2012        
2013           - M{S{xi}(P)=P/3} where M{P<3}           - M{S{xi}(P)=P/3} where M{P<3}
2014           - M{S{xi}(P)=1/2} otherwise           - M{S{xi}(P)=1/2} otherwise
2015    
2016       @param P: Preclet number       @param P: Preclet number
2017       @type P: L{Scalar<escript.Scalar>}       @type P: L{Scalar<escript.Scalar>}
2018       @return: up-wind weightimg factor       @return: up-wind weightimg factor
2019       @rtype: L{Scalar<escript.Scalar>}       @rtype: L{Scalar<escript.Scalar>}
2020       """       """
2021       c=(P-3.).whereNegative()       c=util.whereNegative(P-3.)
2022       return P/6.*c+1./2.*(1.-c)       return P/6.*c+1./2.*(1.-c)
2023    
2024     def HALF(self,P):     def HALF(self,P):
2025       """       """
2026       Predefined function to set value M{1/2} for M{S{xi}}       Predefined function to set value M{1/2} for M{S{xi}}
2027    
2028       @param P: Preclet number       @param P: Preclet number
2029       @type P: L{Scalar<escript.Scalar>}       @type P: L{Scalar<escript.Scalar>}
2030       @return: up-wind weightimg factor       @return: up-wind weightimg factor
2031       @rtype: L{Scalar<escript.Scalar>}       @rtype: L{Scalar<escript.Scalar>}
2032       """       """
2033       return escript.Scalar(0.5,P.getFunctionSpace())       return escript.Scalar(0.5,P.getFunctionSpace())
2034    
2035     def __calculateXi(self,peclet_factor,Z,h):     def __calculateXi(self,peclet_factor,flux,h):
2036         Z_max=util.Lsup(Z)         flux=util.Lsup(flux)
2037         if Z_max>0.:         if flux_max>0.:
2038            return h*self.__xi(Z*peclet_factor)/(Z+Z_max*self.__TOL)            return h*self.__xi(flux*peclet_factor)/(flux+flux_max*self.__TOL)
2039         else:         else:
2040            return 0.            return 0.
2041    
# Line 2034  class AdvectivePDE(LinearPDE): Line 2048  class AdvectivePDE(LinearPDE):
2048           self.__Xi=escript.Scalar(0.,self.getFunctionSpaceForCoefficient("A"))           self.__Xi=escript.Scalar(0.,self.getFunctionSpaceForCoefficient("A"))
2049           if not C.isEmpty() or not B.isEmpty():           if not C.isEmpty() or not B.isEmpty():
2050              if not C.isEmpty() and not B.isEmpty():              if not C.isEmpty() and not B.isEmpty():
2051                  Z2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))                  flux2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))
2052                  if self.getNumEquations()>1:                  if self.getNumEquations()>1:
2053                     if self.getNumSolutions()>1:                     if self.getNumSolutions()>1:
2054                        for i in range(self.getNumEquations()):                        for i in range(self.getNumEquations()):
2055                           for k in range(self.getNumSolutions()):                           for k in range(self.getNumSolutions()):
2056                              for l in range(self.getDim()): Z2+=(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
2057                          # flux=C-util.reorderComponents(B,[0,2,1])
2058                     else:                     else:
2059                        for i in range(self.getNumEquations()):                        for i in range(self.getNumEquations()):
2060                           for l in range(self.getDim()): Z2+=(C[i,l]-B[i,l])**2                           for l in range(self.getDim()): flux2+=(C[i,l]-B[i,l])**2
2061                          # flux=C-B
2062                  else:                  else:
2063                     if self.getNumSolutions()>1:                     if self.getNumSolutions()>1:
2064                        for k in range(self.getNumSolutions()):                        for k in range(self.getNumSolutions()):
2065                           for l in range(self.getDim()): Z2+=(C[k,l]-B[l,k])**2                           for l in range(self.getDim()): flux2+=(C[k,l]-B[l,k])**2
2066                          # flux=C-util.reorderComponents(B,[1,0])
2067                     else:                     else:
2068                        for l in range(self.getDim()): Z2+=(C[l]-B[l])**2                        for l in range(self.getDim()): flux2+=(C[l]-B[l])**2
2069                  length_of_Z=util.sqrt(Z2)                        #flux=C-B
2070                    length_of_flux=util.sqrt(flux2)
2071              elif C.isEmpty():              elif C.isEmpty():
2072                length_of_Z=util.length(B)                length_of_flux=util.length(B)
2073                  #flux=B
2074              else:              else:
2075                length_of_Z=util.length(C)                length_of_flux=util.length(C)
2076                  #flux=C
2077    
2078              Z_max=util.Lsup(length_of_Z)              #length_of_flux=util.length(flux)
2079              if Z_max>0.:              flux_max=util.Lsup(length_of_flux)
2080                if flux_max>0.:
2081                   # length_of_A=util.inner(flux,util.tensormutiply(A,flux))
2082                 length_of_A=util.length(A)                 length_of_A=util.length(A)
2083                 A_max=util.Lsup(length_of_A)                 A_max=util.Lsup(length_of_A)
2084                 if A_max>0:                 if A_max>0:
2085                      inv_A=1./(length_of_A+A_max*self.__TOL)                      inv_A=1./(length_of_A+A_max*self.__TOL)
2086                 else:                 else:
2087                      inv_A=1./self.__TOL                      inv_A=1./self.__TOL
2088                 peclet_number=length_of_Z*h/2*inv_A                 peclet_number=length_of_flux*h/2*inv_A
2089                 xi=self.__xi(peclet_number)                 xi=self.__xi(peclet_number)
2090                 self.__Xi=h*xi/(length_of_Z+Z_max*self.__TOL)                 self.__Xi=h*xi/(length_of_flux+flux_max*self.__TOL)
2091                 self.trace("preclet number = %e"%util.Lsup(peclet_number))                 self.trace("preclet number = %e"%util.Lsup(peclet_number))
2092        return self.__Xi        return self.__Xi
2093    
# Line 2103  class AdvectivePDE(LinearPDE): Line 2125  class AdvectivePDE(LinearPDE):
2125                        for k in range(self.getNumSolutions()):                        for k in range(self.getNumSolutions()):
2126                           for l in range(self.getDim()):                           for l in range(self.getDim()):
2127                              if not C.isEmpty() and not B.isEmpty():                              if not C.isEmpty() and not B.isEmpty():
2128                                   # tmp=C-util.reorderComponents(B,[0,2,1])
2129                                   # Aout=Aout+Xi*util.generalTensorProduct(util.reorder(tmp,[1,2,0]),tmp,offset=1)
2130                                 for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*(C[p,i,j]-B[p,j,i])*(C[p,k,l]-B[p,l,k])                                 for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*(C[p,i,j]-B[p,j,i])*(C[p,k,l]-B[p,l,k])
2131                              elif C.isEmpty():                              elif C.isEmpty():
2132                                 for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*B[p,j,i]*B[p,l,k]                                 for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*B[p,j,i]*B[p,l,k]
2133                                   # Aout=Aout+Xi*util.generalTensorProduct(util.reorder(B,[2,1,0]),util.reorder(B,[0,2,1]),offset=1)
2134                              else:                              else:
2135                                 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]
2136                                   # Aout=Aout+Xi*util.generalTensorProduct(util.reorder(C,[1,2,0]),C,offset=1)
2137              else:              else:
2138                  for j in range(self.getDim()):                  for j in range(self.getDim()):
2139                     for l in range(self.getDim()):                     for l in range(self.getDim()):
# Line 2117  class AdvectivePDE(LinearPDE): Line 2143  class AdvectivePDE(LinearPDE):
2143                            Aout[j,l]+=Xi*B[j]*B[l]                            Aout[j,l]+=Xi*B[j]*B[l]
2144                        else:                        else:
2145                            Aout[j,l]+=Xi*C[j]*C[l]                            Aout[j,l]+=Xi*C[j]*C[l]
2146                     # if not C.isEmpty() and not B.isEmpty():
2147                     #    tmp=C-B
2148                     #    Aout=Aout+Xi*util.outer(tmp,tmp)
2149                     # elif C.isEmpty():
2150                     #    Aout=Aout+Xi*util.outer(B,B)
2151                     # else:
2152                     # Aout=Aout+Xi*util.outer(C,C)
2153           return Aout           return Aout
2154       elif name == "B" :       elif name == "B" :
2155           B=self.getCoefficient("B")           B=self.getCoefficient("B")
# Line 2137  class AdvectivePDE(LinearPDE): Line 2170  class AdvectivePDE(LinearPDE):
2170                       for i in range(self.getNumEquations()):                       for i in range(self.getNumEquations()):
2171                          for j in range(self.getDim()):                          for j in range(self.getDim()):
2172                             Bout[i,j,k]+=tmp*C[p,i,j]                             Bout[i,j,k]+=tmp*C[p,i,j]
2173                               # Bout=Bout+Xi*util.generalTensorProduct(util.reorder(C,[1,2,0]),D,offset=1)
2174              else:              else:
2175                 tmp=Xi*D                 tmp=Xi*D
2176                 for j in range(self.getDim()): Bout[j]+=tmp*C[j]                 for j in range(self.getDim()): Bout[j]+=tmp*C[j]
2177                   # Bout=Bout+Xi*D*C
2178           return Bout           return Bout
2179       elif name == "C" :       elif name == "C" :
2180           B=self.getCoefficient("B")           B=self.getCoefficient("B")
# Line 2160  class AdvectivePDE(LinearPDE): Line 2195  class AdvectivePDE(LinearPDE):
2195                        for i in range(self.getNumEquations()):                        for i in range(self.getNumEquations()):
2196                          for l in range(self.getDim()):                          for l in range(self.getDim()):
2197                                   Cout[i,k,l]+=tmp*B[p,l,i]                                   Cout[i,k,l]+=tmp*B[p,l,i]
2198                                     # Cout=Cout+Xi*B[p,l,i]*D[p,k]
2199              else:              else:
2200                 tmp=Xi*D                 tmp=Xi*D
2201                 for j in range(self.getDim()): Cout[j]+=tmp*B[j]                 for j in range(self.getDim()): Cout[j]+=tmp*B[j]
2202                   # Cout=Cout+tmp*D*B
2203           return Cout           return Cout
2204       elif name == "D" :       elif name == "D" :
2205           return self.getCoefficient("D")           return self.getCoefficient("D")
# Line 2186  class AdvectivePDE(LinearPDE): Line 2223  class AdvectivePDE(LinearPDE):
2223                         for j in range(self.getDim()):                         for j in range(self.getDim()):
2224                            if not C.isEmpty() and not B.isEmpty():                            if not C.isEmpty() and not B.isEmpty():
2225                               Xout[i,j]+=tmp*(C[p,i,j]-B[p,j,i])                               Xout[i,j]+=tmp*(C[p,i,j]-B[p,j,i])
2226                                 # Xout=X_out+Xi*util.inner(Y,C-util.reorderComponents(B,[0,2,1]),offset=1)
2227                            elif C.isEmpty():                            elif C.isEmpty():
2228                               Xout[i,j]-=tmp*B[p,j,i]                               Xout[i,j]-=tmp*B[p,j,i]
2229                                 # Xout=X_out-Xi*util.inner(Y,util.reorderComponents(B,[0,2,1]),offset=1)
2230                            else:                            else:
2231                               Xout[i,j]+=tmp*C[p,i,j]                               Xout[i,j]+=tmp*C[p,i,j]
2232                                 # Xout=X_out+Xi*util.inner(Y,C,offset=1)
2233              else:              else:
2234                   tmp=Xi*Y                   tmp=Xi*Y
2235                   for j in range(self.getDim()):                   for j in range(self.getDim()):
2236                      if not C.isEmpty() and not B.isEmpty():                      if not C.isEmpty() and not B.isEmpty():
2237                         Xout[j]+=tmp*(C[j]-B[j])                         Xout[j]+=tmp*(C[j]-B[j])
2238                           # Xout=Xout+Xi*Y*(C-B)
2239                      elif C.isEmpty():                      elif C.isEmpty():
2240                         Xout[j]-=tmp*B[j]                         Xout[j]-=tmp*B[j]
2241                           # Xout=Xout-Xi*Y*B
2242                      else:                      else:
2243                         Xout[j]+=tmp*C[j]                         Xout[j]+=tmp*C[j]
2244                           # Xout=Xout+Xi*Y*C
2245           return Xout           return Xout
2246       elif name == "Y" :       elif name == "Y" :
2247           return self.getCoefficient("Y")           return self.getCoefficient("Y")
# Line 2217  class AdvectivePDE(LinearPDE): Line 2260  class AdvectivePDE(LinearPDE):
2260       else:       else:
2261          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
2262    
2263    class AdvectionDiffusion(LinearPDE):
2264       """
2265       Class to define PDE equation of the unisotropic advection-diffusion problem, which is genear L{LinearPDE} of the form
2266    
2267       M{S{omega}*u + inner(v,grad(u))- grad(matrixmult(k_bar,grad(u))[j])[j] = f}
2268    
2269       with natural boundary conditons
2270    
2271       M{n[j]*matrixmult(k,grad(u))[j] = g- S{alpha}u }
2272    
2273       and constraints:
2274    
2275       M{u=r} where M{q>0}
2276    
2277       and
2278    
2279       M{k_bar[i,j]=k[i,j]+upwind[i]*upwind[j]}
2280    
2281       """
2282    
2283       def __init__(self,domain,debug=False):
2284         """
2285         initializes a new Poisson equation
2286    
2287         @param domain: domain of the PDE
2288         @type domain: L{Domain<escript.Domain>}
2289         @param debug: if True debug informations are printed.
2290    
2291         """
2292         super(AdvectionDiffusion, self).__init__(domain,1,1,debug)
2293         self.COEFFICIENTS={"omega": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
2294                            "k": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_DIM,PDECoefficient.BY_DIM),PDECoefficient.OPERATOR),
2295                            "f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
2296                            "v": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_DIM,),PDECoefficient.OPERATOR),
2297                            "upwind": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_DIM,),PDECoefficient.OPERATOR),
2298                            "alpha": PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
2299                            "g": PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
2300                            "r": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH),
2301                            "q": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}
2302    
2303       def setValue(self,**coefficients):
2304         """
2305         sets new values to coefficients
2306    
2307         @param coefficients: new values assigned to coefficients
2308         @keyword omega: value for coefficient M{S{omega}}
2309         @type omega: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
2310         @keyword k: value for coefficient M{k}
2311         @type k: any type that can be casted to L{Tensor<escript.Tensor>} object on L{Function<escript.Function>}.
2312         @keyword v: value for coefficient M{v}
2313         @type v: any type that can be casted to L{Vector<escript.Vector>} object on L{Function<escript.Function>}.
2314         @keyword upwind: value for upwind term M{upwind}
2315         @type upwind: any type that can be casted to L{Vector<escript.Vector>} object on L{Function<escript.Function>}.
2316         @keyword f: value for right hand side M{f}
2317         @type f: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
2318         @keyword alpha: value for right hand side M{S{alpha}}
2319         @type alpha: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
2320         @keyword g: value for right hand side M{g}
2321         @type g: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
2322         @keyword r: prescribed values M{r} for the solution in constraints.
2323         @type r: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2324                   depending of reduced order is used for the representation of the equation.
2325         @keyword q: mask for location of constraints
2326         @type q: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2327                   depending of reduced order is used for the representation of the equation.
2328         @raise IllegalCoefficient: if an unknown coefficient keyword is used.
2329         """
2330         super(AdvectionDiffusion, self).setValue(**coefficients)
2331    
2332       def getCoefficientOfGeneralPDE(self,name):
2333         """
2334         return the value of the coefficient name of the general PDE
2335    
2336         @param name: name of the coefficient requested.
2337         @type name: C{string}
2338         @return: the value of the coefficient  name
2339         @rtype: L{Data<escript.Data>}
2340         @raise IllegalCoefficient: if name is not one of coefficients
2341                      "A", M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact}, M{r} or M{q}.
2342         @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
2343         """
2344         if name == "A" :
2345             return self.getCoefficient("k")+util.outer(self.getCoefficient("upwind"),self.getCoefficient("upwind"))
2346         elif name == "B" :
2347             return escript.Data()
2348         elif name == "C" :
2349             return self.getCoefficient("v")
2350         elif name == "D" :
2351             return self.getCoefficient("omega")
2352         elif name == "X" :
2353             return escript.Data()
2354         elif name == "Y" :
2355             return self.getCoefficient("f")
2356         elif name == "d" :
2357             return self.getCoefficient("alpha")
2358         elif name == "y" :
2359             return self.getCoefficient("g")
2360         elif name == "d_contact" :
2361             return escript.Data()
2362         elif name == "y_contact" :
2363             return escript.Data()
2364         elif name == "r" :
2365             return self.getCoefficient("r")
2366         elif name == "q" :
2367             return self.getCoefficient("q")
2368         else:
2369            raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
2370    
2371    
2372  # $Log$  # $Log$
2373    # Revision 1.14  2005/09/22 01:54:57  jgs
2374    # Merge of development branch dev-02 back to main trunk on 2005-09-22
2375    #
2376  # Revision 1.13  2005/09/15 03:44:19  jgs  # Revision 1.13  2005/09/15 03:44:19  jgs
2377  # 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
2378  #  #
# Line 2231  class AdvectivePDE(LinearPDE): Line 2385  class AdvectivePDE(LinearPDE):
2385  # Revision 1.10  2005/08/12 01:45:36  jgs  # Revision 1.10  2005/08/12 01:45:36  jgs
2386  # 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
2387  #  #
2388    # Revision 1.9.2.17  2005/09/21 07:03:33  matt
2389    # PDECoefficient and LinearPDE are now new style classes (introduced in Python
2390    # 2.2). Classes Poisson, Helmholtz, LameEquation and AdvectivePDE have been
2391    # modified to instead use portable/cooperative "super" calls to extend base
2392    # class methods.
2393    #
2394    # Revision 1.9.2.16  2005/09/16 01:54:37  matt
2395    # Removed redundant if-loop.
2396    #
2397  # Revision 1.9.2.15  2005/09/14 08:09:18  matt  # Revision 1.9.2.15  2005/09/14 08:09:18  matt
2398  # Added "REDUCED" solution PDECoefficient descriptors for LinearPDEs.  # Added "REDUCED" solution PDECoefficient descriptors for LinearPDEs.
2399  #  #

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

  ViewVC Help
Powered by ViewVC 1.1.26