/[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

trunk/esys2/escript/py_src/linearPDEs.py revision 149 by jgs, Thu Sep 1 03:31:39 2005 UTC trunk/escript/py_src/linearPDEs.py revision 425 by gross, Tue Jan 10 04:10:39 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 62  class PDECoefficient: Line 62  class PDECoefficient:
62      @cvar BOUNDARY: indicator that coefficient is defined on the boundary of the PDE domain      @cvar BOUNDARY: indicator that coefficient is defined on the boundary of the PDE domain
63      @cvar CONTACT: indicator that coefficient is defined on the contact region within the PDE domain      @cvar CONTACT: indicator that coefficient is defined on the contact region within the PDE domain
64      @cvar SOLUTION: indicator that coefficient is defined trough a solution of the PDE      @cvar SOLUTION: indicator that coefficient is defined trough a solution of the PDE
65        @cvar REDUCED: indicator that coefficient is defined trough a reduced solution of the PDE
66      @cvar BY_EQUATION: indicator that the dimension of the coefficient shape is defined by the number PDE equations      @cvar BY_EQUATION: indicator that the dimension of the coefficient shape is defined by the number PDE equations
67      @cvar BY_SOLUTION: indicator that the dimension of the coefficient shape is defined by the number PDE solutions      @cvar BY_SOLUTION: indicator that the dimension of the coefficient shape is defined by the number PDE solutions
68      @cvar BY_DIM: indicator that the dimension of the coefficient shape is defined by the spatial dimension      @cvar BY_DIM: indicator that the dimension of the coefficient shape is defined by the spatial dimension
# Line 74  class PDECoefficient: Line 75  class PDECoefficient:
75      BOUNDARY=1      BOUNDARY=1
76      CONTACT=2      CONTACT=2
77      SOLUTION=3      SOLUTION=3
78        REDUCED=4
79      BY_EQUATION=5      BY_EQUATION=5
80      BY_SOLUTION=6      BY_SOLUTION=6
81      BY_DIM=7      BY_DIM=7
# Line 84  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}         @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 99  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 123  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:
141                return escript.ReducedSolution(domain)
142    
143      def getValue(self):      def getValue(self):
144         """         """
# Line 306  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 336  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 356  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 365  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 380  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 392  class LinearPDE: Line 397  class LinearPDE:
397     In this case the the coefficient M{d_contact} and M{y_contact} are eaach scalar     In this case the the coefficient M{d_contact} and M{y_contact} are eaach scalar
398     both in the L{FunctionOnContactZero<escript.FunctionOnContactZero>} or L{FunctionOnContactOne<escript.FunctionOnContactOne>}.     both in the L{FunctionOnContactZero<escript.FunctionOnContactZero>} or L{FunctionOnContactOne<escript.FunctionOnContactOne>}.
399    
400       @cvar DEFAULT: The default method used to solve the system of linear equations
    @cvar DEFAULT_METHOD: The default method used to solve the system of linear equations  
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 410  class LinearPDE: Line 414  class LinearPDE:
414     @cvar NO_REORDERING: No matrix reordering allowed     @cvar NO_REORDERING: No matrix reordering allowed
415     @cvar MINIMUM_FILL_IN: Reorder matrix to reduce fill-in during factorization     @cvar MINIMUM_FILL_IN: Reorder matrix to reduce fill-in during factorization
416     @cvar NESTED_DISSECTION: Reorder matrix to improve load balancing during factorization     @cvar NESTED_DISSECTION: Reorder matrix to improve load balancing during factorization
417       @cvar PASO: PASO solver package
418       @cvar SCSL: SGI SCSL solver library
419       @cvar MKL: Intel's MKL solver library
420       @cvar UMFPACK: the UMFPACK library
421       @cvar ITERATIVE: The default iterative solver
422    
423     """     """
424     DEFAULT_METHOD=0     DEFAULT= 0
425     DIRECT=1     DIRECT= 1
426     CHOLEVSKY=2     CHOLEVSKY= 2
427     PCG=3     PCG= 3
428     CR=4     CR= 4
429     CGS=5     CGS= 5
430     BICGSTAB=6     BICGSTAB= 6
431     SSOR=7     SSOR= 7
432     ILU0=8     ILU0= 8
433     ILUT=9     ILUT= 9
434     JACOBI=10     JACOBI= 10
435     GMRES=11     GMRES= 11
436     PRES20=12     PRES20= 12
437     LUMPING=13     LUMPING= 13
438     NO_REORDERING=0     NO_REORDERING= 17
439     MINIMUM_FILL_IN=1     MINIMUM_FILL_IN= 18
440     NESTED_DISSECTION=2     NESTED_DISSECTION= 19
441       SCSL= 14
442       MKL= 15
443       UMFPACK= 16
444       ITERATIVE= 20
445       PASO= 21
446    
447     __TOL=1.e-13     __TOL=1.e-13
448       __PACKAGE_KEY="package"
449     __METHOD_KEY="method"     __METHOD_KEY="method"
450     __SYMMETRY_KEY="symmetric"     __SYMMETRY_KEY="symmetric"
451     __TOLERANCE_KEY="tolerance"     __TOLERANCE_KEY="tolerance"
452       __PRECONDITIONER_KEY="preconditioner"
453    
454    
455     def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):     def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):
# Line 448  class LinearPDE: Line 465  class LinearPDE:
465       @param debug: if True debug informations are printed.       @param debug: if True debug informations are printed.
466    
467       """       """
468         super(LinearPDE, self).__init__()
469       #       #
470       #   the coefficients of the general PDE:       #   the coefficients of the general PDE:
471       #       #
# Line 479  class LinearPDE: Line 497  class LinearPDE:
497       self.__reduce_equation_order=False       self.__reduce_equation_order=False
498       self.__reduce_solution_order=False       self.__reduce_solution_order=False
499       self.__tolerance=1.e-8       self.__tolerance=1.e-8
500       self.__solver_method=self.DEFAULT_METHOD       self.__solver_method=self.DEFAULT
501       self.__matrix_type=self.__domain.getSystemMatrixTypeId(self.DEFAULT_METHOD,False)       self.__solver_package=self.DEFAULT
502         self.__preconditioner=self.DEFAULT
503         self.__matrix_type=self.__domain.getSystemMatrixTypeId(self.DEFAULT,self.DEFAULT,False)
504       self.__sym=False       self.__sym=False
505    
506       self.resetCoefficients()       self.resetCoefficients()
# Line 583  class LinearPDE: Line 603  class LinearPDE:
603       @rtype: L{bool}       @rtype: L{bool}
604       """       """
605       return self.__reduce_solution_order       return self.__reduce_solution_order
606    
607     def getFunctionSpaceForEquation(self):     def getFunctionSpaceForEquation(self):
608       """       """
609       returns the L{FunctionSpace<escript.FunctionSpace>} used to discretize the equation       returns the L{FunctionSpace<escript.FunctionSpace>} used to discretize the equation
# Line 754  class LinearPDE: Line 774  class LinearPDE:
774         @type verbose: C{bool}         @type verbose: C{bool}
775         @keyword reordering: reordering scheme to be used during elimination. Allowed values are         @keyword reordering: reordering scheme to be used during elimination. Allowed values are
776                              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}  
777         @keyword iter_max: maximum number of iteration steps allowed.         @keyword iter_max: maximum number of iteration steps allowed.
778         @keyword drop_tolerance: threshold for drupping in L{ILUT}         @keyword drop_tolerance: threshold for drupping in L{ILUT}
779         @keyword drop_storage: maximum of allowed memory in L{ILUT}         @keyword drop_storage: maximum of allowed memory in L{ILUT}
# Line 768  class LinearPDE: Line 786  class LinearPDE:
786               self.__solution=self.copyConstraint(f*mat)               self.__solution=self.copyConstraint(f*mat)
787            else:            else:
788               options[self.__TOLERANCE_KEY]=self.getTolerance()               options[self.__TOLERANCE_KEY]=self.getTolerance()
789               options[self.__METHOD_KEY]=self.getSolverMethod()               options[self.__METHOD_KEY]=self.getSolverMethod()[0]
790                 options[self.__PRECONDITIONER_KEY]=self.getSolverMethod()[1]
791                 options[self.__PACKAGE_KEY]=self.getSolverPackage()
792               options[self.__SYMMETRY_KEY]=self.isSymmetric()               options[self.__SYMMETRY_KEY]=self.isSymmetric()
793               self.trace("PDE is resolved.")               self.trace("PDE is resolved.")
794               self.trace("solver options: %s"%str(options))               self.trace("solver options: %s"%str(options))
# Line 782  class LinearPDE: Line 802  class LinearPDE:
802    
803       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]}
804    
805       or       or
806    
807       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]}
808    
# Line 796  class LinearPDE: Line 816  class LinearPDE:
816     # =============================================================================     # =============================================================================
817     #   solver settings:     #   solver settings:
818     # =============================================================================     # =============================================================================
819     def setSolverMethod(self,solver=None):     def setSolverMethod(self,solver=None,preconditioner=None):
820         """         """
821         sets a new solver         sets a new solver
822    
823         @param solver: sets a new solver method.         @param solver: sets a new solver method.
824         @type solver: one of L{DEFAULT_METHOD}, 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}.
825           @param preconditioner: sets a new solver method.
826         """         @type solver: one of L{DEFAULT}, L{JACOBI} L{ILU0}, L{ILUT},L{SSOR}
827         if solver==None: solve=self.DEFAULT_METHOD         """
828         if not solver==self.getSolverMethod():         if solver==None: solve=self.DEFAULT
829           if preconditioner==None: preconditioner=self.DEFAULT
830           if not (solver,preconditioner)==self.getSolverMethod():
831             self.__solver_method=solver             self.__solver_method=solver
832               self.__preconditioner=preconditioner
833             self.__checkMatrixType()             self.__checkMatrixType()
834             self.trace("New solver is %s"%self.getSolverMethodName())             self.trace("New solver is %s"%self.getSolverMethodName())
835    
# Line 819  class LinearPDE: Line 842  class LinearPDE:
842         """         """
843    
844         m=self.getSolverMethod()         m=self.getSolverMethod()
845         if m==self.DEFAULT_METHOD: return "DEFAULT_METHOD"         p=self.getSolverPackage()
846         elif m==self.DIRECT: return "DIRECT"         method=""
847         elif m==self.CHOLEVSKY: return "CHOLEVSKY"         if m[0]==self.DEFAULT: method="DEFAULT"
848         elif m==self.PCG: return "PCG"         elif m[0]==self.DIRECT: method= "DIRECT"
849         elif m==self.CR: return "CR"         elif m[0]==self.ITERATIVE: method= "ITERATIVE"
850         elif m==self.CGS: return "CGS"         elif m[0]==self.CHOLEVSKY: method= "CHOLEVSKY"
851         elif m==self.BICGSTAB: return "BICGSTAB"         elif m[0]==self.PCG: method= "PCG"
852         elif m==self.SSOR: return "SSOR"         elif m[0]==self.CR: method= "CR"
853         elif m==self.GMRES: return "GMRES"         elif m[0]==self.CGS: method= "CGS"
854         elif m==self.PRES20: return "PRES20"         elif m[0]==self.BICGSTAB: method= "BICGSTAB"
855         elif m==self.LUMPING: return "LUMPING"         elif m[0]==self.SSOR: method= "SSOR"
856         return None         elif m[0]==self.GMRES: method= "GMRES"
857           elif m[0]==self.PRES20: method= "PRES20"
858           elif m[0]==self.LUMPING: method= "LUMPING"
859           if m[1]==self.DEFAULT: method+="+DEFAULT"
860           elif m[1]==self.JACOBI: method+= "+JACOBI"
861           elif m[1]==self.ILU0: method+= "+ILU0"
862           elif m[1]==self.ILUT: method+= "+ILUT"
863           elif m[1]==self.SSOR: method+= "+SSOR"
864           if p==self.DEFAULT: package="DEFAULT"
865           elif p==self.PASO: package= "PASO"
866           elif p==self.MKL: package= "MKL"
867           elif p==self.SCSL: package= "SCSL"
868           elif p==self.UMFPACK: package= "UMFPACK"
869           else : method="unknown"
870           return "%s solver of %s package"%(method,package)
871    
872    
873     def getSolverMethod(self):     def getSolverMethod(self):
874         """         """
875         returns the solver method         returns the solver method
876    
877         @return: the solver method currently be used.         @return: the solver method currently be used.
878         @rtype: C{int}         @rtype: C{int}
879         """         """
880         return self.__solver_method         return self.__solver_method,self.__preconditioner
881    
882       def setSolverPackage(self,package=None):
883           """
884           sets a new solver package
885    
886           @param solver: sets a new solver method.
887           @type solver: one of L{DEFAULT}, L{PASO} L{SCSL}, L{MKL}, L{UMLPACK}
888           """
889           if package==None: package=self.DEFAULT
890           if not package==self.getSolverPackage():
891               self.__solver_method=solver
892               self.__checkMatrixType()
893               self.trace("New solver is %s"%self.getSolverMethodName())
894    
895       def getSolverPackage(self):
896           """
897           returns the package of the solver
898    
899           @return: the solver package currently being used.
900           @rtype: C{int}
901           """
902           return self.__solver_package
903    
904     def isUsingLumping(self):     def isUsingLumping(self):
905        """        """
# Line 849  class LinearPDE: Line 908  class LinearPDE:
908        @return: True is lumping is currently used a solver method.        @return: True is lumping is currently used a solver method.
909        @rtype: C{bool}        @rtype: C{bool}
910        """        """
911        return self.getSolverMethod()==self.LUMPING        return self.getSolverMethod()[0]==self.LUMPING
912    
913     def setTolerance(self,tol=1.e-8):     def setTolerance(self,tol=1.e-8):
914         """         """
# Line 1042  class LinearPDE: Line 1101  class LinearPDE:
1101       """       """
1102       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.
1103       """       """
1104       new_matrix_type=self.getDomain().getSystemMatrixTypeId(self.getSolverMethod(),self.isSymmetric())       new_matrix_type=self.getDomain().getSystemMatrixTypeId(self.getSolverMethod()[0],self.getSolverPackage(),self.isSymmetric())
1105       if not new_matrix_type==self.__matrix_type:       if not new_matrix_type==self.__matrix_type:
1106           self.trace("Matrix type is now %d."%new_matrix_type)           self.trace("Matrix type is now %d."%new_matrix_type)
1107           self.__matrix_type=new_matrix_type           self.__matrix_type=new_matrix_type
# Line 1555  class Poisson(LinearPDE): Line 1614  class Poisson(LinearPDE):
1614    
1615     """     """
1616    
1617     def __init__(self,domain,f=escript.Data(),q=escript.Data(),debug=False):     def __init__(self,domain,debug=False):
1618       """       """
1619       initializes a new Poisson equation       initializes a new Poisson equation
1620    
# Line 1564  class Poisson(LinearPDE): Line 1623  class Poisson(LinearPDE):
1623       @param debug: if True debug informations are printed.       @param debug: if True debug informations are printed.
1624    
1625       """       """
1626       LinearPDE.__init__(self,domain,1,1,debug)       super(Poisson, self).__init__(domain,1,1,debug)
1627       self.COEFFICIENTS={"f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),       self.COEFFICIENTS={"f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1628                            "q": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}                            "q": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}
1629       self.setSymmetryOn()       self.setSymmetryOn()
# Line 1581  class Poisson(LinearPDE): Line 1640  class Poisson(LinearPDE):
1640                 depending of reduced order is used for the representation of the equation.                 depending of reduced order is used for the representation of the equation.
1641       @raise IllegalCoefficient: if an unknown coefficient keyword is used.       @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1642       """       """
1643       LinearPDE.setValue(self,**coefficients)       super(Poisson, self).setValue(**coefficients)
1644    
1645     def getCoefficientOfGeneralPDE(self,name):     def getCoefficientOfGeneralPDE(self,name):
1646       """       """
# Line 1595  class Poisson(LinearPDE): Line 1654  class Poisson(LinearPDE):
1654       @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.
1655       """       """
1656       if name == "A" :       if name == "A" :
1657           return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))           return escript.Data(util.kronecker(self.getDim()),escript.Function(self.getDomain()))
1658       elif name == "B" :       elif name == "B" :
1659           return escript.Data()           return escript.Data()
1660       elif name == "C" :       elif name == "C" :
# Line 1646  class Helmholtz(LinearPDE): Line 1705  class Helmholtz(LinearPDE):
1705       @param debug: if True debug informations are printed.       @param debug: if True debug informations are printed.
1706    
1707       """       """
1708       LinearPDE.__init__(self,domain,1,1,debug)       super(Helmholtz, self).__init__(domain,1,1,debug)
1709       self.COEFFICIENTS={"omega": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),       self.COEFFICIENTS={"omega": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
1710                          "k": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),                          "k": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
1711                          "f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),                          "f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
# Line 1679  class Helmholtz(LinearPDE): Line 1738  class Helmholtz(LinearPDE):
1738                 depending of reduced order is used for the representation of the equation.                 depending of reduced order is used for the representation of the equation.
1739       @raise IllegalCoefficient: if an unknown coefficient keyword is used.       @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1740       """       """
1741       LinearPDE.setValue(self,**coefficients)       super(Helmholtz, self).setValue(**coefficients)
1742    
1743     def getCoefficientOfGeneralPDE(self,name):     def getCoefficientOfGeneralPDE(self,name):
1744       """       """
# Line 1737  class LameEquation(LinearPDE): Line 1796  class LameEquation(LinearPDE):
1796     """     """
1797    
1798     def __init__(self,domain,debug=False):     def __init__(self,domain,debug=False):
1799         LinearPDE.__init__(self,domain,domain.getDim(),domain.getDim(),debug)        super(LameEquation, self).__init__(domain,\
1800         self.COEFFICIENTS={ "lame_lambda"  : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),                                           domain.getDim(),domain.getDim(),debug)
1801          self.COEFFICIENTS={ "lame_lambda"  : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),
1802                            "lame_mu"      : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),                            "lame_mu"      : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),
1803                            "F"            : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),                            "F"            : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1804                            "sigma"        : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM),PDECoefficient.RIGHTHANDSIDE),                            "sigma"        : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM),PDECoefficient.RIGHTHANDSIDE),
1805                            "f"            : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),                            "f"            : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1806                            "r"            : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH),                            "r"            : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH),
1807                            "q"            : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}                            "q"            : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}
1808         self.setSymmetryOn()        self.setSymmetryOn()
1809    
1810     def setValue(self,**coefficients):     def setValue(self,**coefficients):
1811       """       """
# Line 1770  class LameEquation(LinearPDE): Line 1830  class LameEquation(LinearPDE):
1830                 depending of reduced order is used for the representation of the equation.                 depending of reduced order is used for the representation of the equation.
1831       @raise IllegalCoefficient: if an unknown coefficient keyword is used.       @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1832       """       """
1833       LinearPDE.setValue(self,**coefficients)       super(LameEquation, self).setValue(**coefficients)
1834    
1835     def getCoefficientOfGeneralPDE(self,name):     def getCoefficientOfGeneralPDE(self,name):
1836       """       """
# Line 1826  class AdvectivePDE(LinearPDE): Line 1886  class AdvectivePDE(LinearPDE):
1886    
1887     M{Z[j]=C[j]-B[j]}     M{Z[j]=C[j]-B[j]}
1888    
1889     or     or
1890    
1891     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]}
1892    
1893     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
1894     X{Pelclet number} M{P} is used. It is defined as     X{Pelclet number} M{P} is used. It is defined as
1895    
1896     M{P=h|Z|/(2|A|)}     M{P=h|Z|/(2|A|)}
1897    
1898     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
1899     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}}.
1900    
1901     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 1863  class AdvectivePDE(LinearPDE): Line 1923  class AdvectivePDE(LinearPDE):
1923     """     """
1924     def __init__(self,domain,numEquations=None,numSolutions=None,xi=None,debug=False):     def __init__(self,domain,numEquations=None,numSolutions=None,xi=None,debug=False):
1925        """        """
1926        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>}
1927    
1928        @param domain: domain of the PDE        @param domain: domain of the PDE
1929        @type domain: L{Domain<escript.Domain>}        @type domain: L{Domain<escript.Domain>}
# Line 1871  class AdvectivePDE(LinearPDE): Line 1931  class AdvectivePDE(LinearPDE):
1931                             is exracted from the PDE coefficients.                             is exracted from the PDE coefficients.
1932        @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
1933                             is exracted from the PDE coefficients.                             is exracted from the PDE coefficients.
1934        @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
1935                   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.
1936        @type xi: callable object which returns a L{Scalar<escript.Scalar>} object.        @type xi: callable object which returns a L{Scalar<escript.Scalar>} object.
1937        @param debug: if True debug informations are printed.        @param debug: if True debug informations are printed.
1938        """        """
1939          super(AdvectivePDE, self).__init__(domain,\
1940        LinearPDE.__init__(self,domain,numEquations,numSolutions,debug)                                           numEquations,numSolutions,debug)
1941        if xi==None:        if xi==None:
1942           self.__xi=AdvectivePDE.ELMAN_RAMAGE           self.__xi=AdvectivePDE.ELMAN_RAMAGE
1943        else:        else:
1944           self.__xi=xi           self.__xi=xi
1945        self.__Xi=escript.Data()        self.__Xi=escript.Data()
1946    
1947     def setValue(self,**coefficients):     def setValue(**coefficients):
1948        """        """
1949        sets new values to coefficients        sets new values to coefficients
1950    
# Line 1921  class AdvectivePDE(LinearPDE): Line 1981  class AdvectivePDE(LinearPDE):
1981    
1982        """        """
1983        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()
1984        LinearPDE.setValue(self,**coefficients)        super(AdvectivePDE, self).setValue(**coefficients)
1985    
1986     def ELMAN_RAMAGE(self,P):     def ELMAN_RAMAGE(self,P):
1987       """       """
1988       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}.
1989       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)
1990            - M{S{xi}(P)=0} for M{P<1}            - M{S{xi}(P)=0} for M{P<1}
1991            - M{S{xi}(P)=(1-1/P)/2} otherwise            - M{S{xi}(P)=(1-1/P)/2} otherwise
1992    
1993       @param P: Preclet number       @param P: Preclet number
1994       @type P: L{Scalar<escript.Scalar>}       @type P: L{Scalar<escript.Scalar>}
1995       @return: up-wind weightimg factor       @return: up-wind weightimg factor
1996       @rtype: L{Scalar<escript.Scalar>}       @rtype: L{Scalar<escript.Scalar>}
1997       """       """
1998       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))
1999    
2000     def SIMPLIFIED_BROOK_HUGHES(self,P):     def SIMPLIFIED_BROOK_HUGHES(self,P):
2001       """       """
2002       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}.
2003       The original methods is       The original methods is
2004        
2005       M{S{xi}(P)=coth(P)-1/P}       M{S{xi}(P)=coth(P)-1/P}
2006    
2007       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:
2008        
2009           - M{S{xi}(P)=P/3} where M{P<3}           - M{S{xi}(P)=P/3} where M{P<3}
2010           - M{S{xi}(P)=1/2} otherwise           - M{S{xi}(P)=1/2} otherwise
2011    
2012       @param P: Preclet number       @param P: Preclet number
2013       @type P: L{Scalar<escript.Scalar>}       @type P: L{Scalar<escript.Scalar>}
2014       @return: up-wind weightimg factor       @return: up-wind weightimg factor
2015       @rtype: L{Scalar<escript.Scalar>}       @rtype: L{Scalar<escript.Scalar>}
2016       """       """
2017       c=(P-3.).whereNegative()       c=util.whereNegative(P-3.)
2018       return P/6.*c+1./2.*(1.-c)       return P/6.*c+1./2.*(1.-c)
2019    
2020     def HALF(self,P):     def HALF(self,P):
2021       """       """
2022       Predefined function to set value M{1/2} for M{S{xi}}       Predefined function to set value M{1/2} for M{S{xi}}
2023    
2024       @param P: Preclet number       @param P: Preclet number
2025       @type P: L{Scalar<escript.Scalar>}       @type P: L{Scalar<escript.Scalar>}
2026       @return: up-wind weightimg factor       @return: up-wind weightimg factor
2027       @rtype: L{Scalar<escript.Scalar>}       @rtype: L{Scalar<escript.Scalar>}
2028       """       """
2029       return escript.Scalar(0.5,P.getFunctionSpace())       return escript.Scalar(0.5,P.getFunctionSpace())
2030    
2031     def __calculateXi(self,peclet_factor,Z,h):     def __calculateXi(self,peclet_factor,flux,h):
2032         Z_max=util.Lsup(Z)         flux=util.Lsup(flux)
2033         if Z_max>0.:         if flux_max>0.:
2034            return h*self.__xi(Z*peclet_factor)/(Z+Z_max*self.__TOL)            return h*self.__xi(flux*peclet_factor)/(flux+flux_max*self.__TOL)
2035         else:         else:
2036            return 0.            return 0.
2037    
# Line 1984  class AdvectivePDE(LinearPDE): Line 2044  class AdvectivePDE(LinearPDE):
2044           self.__Xi=escript.Scalar(0.,self.getFunctionSpaceForCoefficient("A"))           self.__Xi=escript.Scalar(0.,self.getFunctionSpaceForCoefficient("A"))
2045           if not C.isEmpty() or not B.isEmpty():           if not C.isEmpty() or not B.isEmpty():
2046              if not C.isEmpty() and not B.isEmpty():              if not C.isEmpty() and not B.isEmpty():
2047                  Z2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))                  flux2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))
2048                  if self.getNumEquations()>1:                  if self.getNumEquations()>1:
2049                     if self.getNumSolutions()>1:                     if self.getNumSolutions()>1:
2050                        for i in range(self.getNumEquations()):                        for i in range(self.getNumEquations()):
2051                           for k in range(self.getNumSolutions()):                           for k in range(self.getNumSolutions()):
2052                              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
2053                          # flux=C-util.reorderComponents(B,[0,2,1])
2054                     else:                     else:
2055                        for i in range(self.getNumEquations()):                        for i in range(self.getNumEquations()):
2056                           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
2057                          # flux=C-B
2058                  else:                  else:
2059                     if self.getNumSolutions()>1:                     if self.getNumSolutions()>1:
2060                        for k in range(self.getNumSolutions()):                        for k in range(self.getNumSolutions()):
2061                           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
2062                          # flux=C-util.reorderComponents(B,[1,0])
2063                     else:                     else:
2064                        for l in range(self.getDim()): Z2+=(C[l]-B[l])**2                        for l in range(self.getDim()): flux2+=(C[l]-B[l])**2
2065                  length_of_Z=util.sqrt(Z2)                        #flux=C-B
2066                    length_of_flux=util.sqrt(flux2)
2067              elif C.isEmpty():              elif C.isEmpty():
2068                length_of_Z=util.length(B)                length_of_flux=util.length(B)
2069                  #flux=B
2070              else:              else:
2071                length_of_Z=util.length(C)                length_of_flux=util.length(C)
2072                  #flux=C
2073    
2074              Z_max=util.Lsup(length_of_Z)              #length_of_flux=util.length(flux)
2075              if Z_max>0.:              flux_max=util.Lsup(length_of_flux)
2076                if flux_max>0.:
2077                   # length_of_A=util.inner(flux,util.tensormutiply(A,flux))
2078                 length_of_A=util.length(A)                 length_of_A=util.length(A)
2079                 A_max=util.Lsup(length_of_A)                 A_max=util.Lsup(length_of_A)
2080                 if A_max>0:                 if A_max>0:
2081                      inv_A=1./(length_of_A+A_max*self.__TOL)                      inv_A=1./(length_of_A+A_max*self.__TOL)
2082                 else:                 else:
2083                      inv_A=1./self.__TOL                      inv_A=1./self.__TOL
2084                 peclet_number=length_of_Z*h/2*inv_A                 peclet_number=length_of_flux*h/2*inv_A
2085                 xi=self.__xi(peclet_number)                 xi=self.__xi(peclet_number)
2086                 self.__Xi=h*xi/(length_of_Z+Z_max*self.__TOL)                 self.__Xi=h*xi/(length_of_flux+flux_max*self.__TOL)
2087                 self.trace("preclet number = %e"%util.Lsup(peclet_number))                 self.trace("preclet number = %e"%util.Lsup(peclet_number))
2088        return self.__Xi        return self.__Xi
2089    
# Line 2053  class AdvectivePDE(LinearPDE): Line 2121  class AdvectivePDE(LinearPDE):
2121                        for k in range(self.getNumSolutions()):                        for k in range(self.getNumSolutions()):
2122                           for l in range(self.getDim()):                           for l in range(self.getDim()):
2123                              if not C.isEmpty() and not B.isEmpty():                              if not C.isEmpty() and not B.isEmpty():
2124                                   # tmp=C-util.reorderComponents(B,[0,2,1])
2125                                   # Aout=Aout+Xi*util.generalTensorProduct(util.reorder(tmp,[1,2,0]),tmp,offset=1)
2126                                 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])
2127                              elif C.isEmpty():                              elif C.isEmpty():
2128                                 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]
2129                                   # Aout=Aout+Xi*util.generalTensorProduct(util.reorder(B,[2,1,0]),util.reorder(B,[0,2,1]),offset=1)
2130                              else:                              else:
2131                                 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]
2132                                   # Aout=Aout+Xi*util.generalTensorProduct(util.reorder(C,[1,2,0]),C,offset=1)
2133              else:              else:
2134                  for j in range(self.getDim()):                  for j in range(self.getDim()):
2135                     for l in range(self.getDim()):                     for l in range(self.getDim()):
# Line 2067  class AdvectivePDE(LinearPDE): Line 2139  class AdvectivePDE(LinearPDE):
2139                            Aout[j,l]+=Xi*B[j]*B[l]                            Aout[j,l]+=Xi*B[j]*B[l]
2140                        else:                        else:
2141                            Aout[j,l]+=Xi*C[j]*C[l]                            Aout[j,l]+=Xi*C[j]*C[l]
2142                     # if not C.isEmpty() and not B.isEmpty():
2143                     #    tmp=C-B
2144                     #    Aout=Aout+Xi*util.outer(tmp,tmp)
2145                     # elif C.isEmpty():
2146                     #    Aout=Aout+Xi*util.outer(B,B)
2147                     # else:
2148                     # Aout=Aout+Xi*util.outer(C,C)
2149           return Aout           return Aout
2150       elif name == "B" :       elif name == "B" :
2151           B=self.getCoefficient("B")           B=self.getCoefficient("B")
# Line 2087  class AdvectivePDE(LinearPDE): Line 2166  class AdvectivePDE(LinearPDE):
2166                       for i in range(self.getNumEquations()):                       for i in range(self.getNumEquations()):
2167                          for j in range(self.getDim()):                          for j in range(self.getDim()):
2168                             Bout[i,j,k]+=tmp*C[p,i,j]                             Bout[i,j,k]+=tmp*C[p,i,j]
2169                               # Bout=Bout+Xi*util.generalTensorProduct(util.reorder(C,[1,2,0]),D,offset=1)
2170              else:              else:
2171                 tmp=Xi*D                 tmp=Xi*D
2172                 for j in range(self.getDim()): Bout[j]+=tmp*C[j]                 for j in range(self.getDim()): Bout[j]+=tmp*C[j]
2173                   # Bout=Bout+Xi*D*C
2174           return Bout           return Bout
2175       elif name == "C" :       elif name == "C" :
2176           B=self.getCoefficient("B")           B=self.getCoefficient("B")
# Line 2110  class AdvectivePDE(LinearPDE): Line 2191  class AdvectivePDE(LinearPDE):
2191                        for i in range(self.getNumEquations()):                        for i in range(self.getNumEquations()):
2192                          for l in range(self.getDim()):                          for l in range(self.getDim()):
2193                                   Cout[i,k,l]+=tmp*B[p,l,i]                                   Cout[i,k,l]+=tmp*B[p,l,i]
2194                                     # Cout=Cout+Xi*B[p,l,i]*D[p,k]
2195              else:              else:
2196                 tmp=Xi*D                 tmp=Xi*D
2197                 for j in range(self.getDim()): Cout[j]+=tmp*B[j]                 for j in range(self.getDim()): Cout[j]+=tmp*B[j]
2198                   # Cout=Cout+tmp*D*B
2199           return Cout           return Cout
2200       elif name == "D" :       elif name == "D" :
2201           return self.getCoefficient("D")           return self.getCoefficient("D")
# Line 2136  class AdvectivePDE(LinearPDE): Line 2219  class AdvectivePDE(LinearPDE):
2219                         for j in range(self.getDim()):                         for j in range(self.getDim()):
2220                            if not C.isEmpty() and not B.isEmpty():                            if not C.isEmpty() and not B.isEmpty():
2221                               Xout[i,j]+=tmp*(C[p,i,j]-B[p,j,i])                               Xout[i,j]+=tmp*(C[p,i,j]-B[p,j,i])
2222                                 # Xout=X_out+Xi*util.inner(Y,C-util.reorderComponents(B,[0,2,1]),offset=1)
2223                            elif C.isEmpty():                            elif C.isEmpty():
2224                               Xout[i,j]-=tmp*B[p,j,i]                               Xout[i,j]-=tmp*B[p,j,i]
2225                                 # Xout=X_out-Xi*util.inner(Y,util.reorderComponents(B,[0,2,1]),offset=1)
2226                            else:                            else:
2227                               Xout[i,j]+=tmp*C[p,i,j]                               Xout[i,j]+=tmp*C[p,i,j]
2228                                 # Xout=X_out+Xi*util.inner(Y,C,offset=1)
2229              else:              else:
2230                   tmp=Xi*Y                   tmp=Xi*Y
2231                   for j in range(self.getDim()):                   for j in range(self.getDim()):
2232                      if not C.isEmpty() and not B.isEmpty():                      if not C.isEmpty() and not B.isEmpty():
2233                         Xout[j]+=tmp*(C[j]-B[j])                         Xout[j]+=tmp*(C[j]-B[j])
2234                           # Xout=Xout+Xi*Y*(C-B)
2235                      elif C.isEmpty():                      elif C.isEmpty():
2236                         Xout[j]-=tmp*B[j]                         Xout[j]-=tmp*B[j]
2237                           # Xout=Xout-Xi*Y*B
2238                      else:                      else:
2239                         Xout[j]+=tmp*C[j]                         Xout[j]+=tmp*C[j]
2240                           # Xout=Xout+Xi*Y*C
2241           return Xout           return Xout
2242       elif name == "Y" :       elif name == "Y" :
2243           return self.getCoefficient("Y")           return self.getCoefficient("Y")
# Line 2167  class AdvectivePDE(LinearPDE): Line 2256  class AdvectivePDE(LinearPDE):
2256       else:       else:
2257          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
2258    
2259    class AdvectionDiffusion(LinearPDE):
2260       """
2261       Class to define PDE equation of the unisotropic advection-diffusion problem, which is genear L{LinearPDE} of the form
2262    
2263       M{S{omega}*u + inner(v,grad(u))- grad(matrixmult(k_bar,grad(u))[j])[j] = f}
2264    
2265       with natural boundary conditons
2266    
2267       M{n[j]*matrixmult(k,grad(u))[j] = g- S{alpha}u }
2268    
2269       and constraints:
2270    
2271       M{u=r} where M{q>0}
2272    
2273       and
2274    
2275       M{k_bar[i,j]=k[i,j]+upwind[i]*upwind[j]}
2276    
2277       """
2278    
2279       def __init__(self,domain,debug=False):
2280         """
2281         initializes a new Poisson equation
2282    
2283         @param domain: domain of the PDE
2284         @type domain: L{Domain<escript.Domain>}
2285         @param debug: if True debug informations are printed.
2286    
2287         """
2288         super(AdvectionDiffusion, self).__init__(domain,1,1,debug)
2289         self.COEFFICIENTS={"omega": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
2290                            "k": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_DIM,PDECoefficient.BY_DIM),PDECoefficient.OPERATOR),
2291                            "f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
2292                            "v": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_DIM,),PDECoefficient.OPERATOR),
2293                            "upwind": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_DIM,),PDECoefficient.OPERATOR),
2294                            "alpha": PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
2295                            "g": PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
2296                            "r": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH),
2297                            "q": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}
2298    
2299       def setValue(self,**coefficients):
2300         """
2301         sets new values to coefficients
2302    
2303         @param coefficients: new values assigned to coefficients
2304         @keyword omega: value for coefficient M{S{omega}}
2305         @type omega: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
2306         @keyword k: value for coefficient M{k}
2307         @type k: any type that can be casted to L{Tensor<escript.Tensor>} object on L{Function<escript.Function>}.
2308         @keyword v: value for coefficient M{v}
2309         @type v: any type that can be casted to L{Vector<escript.Vector>} object on L{Function<escript.Function>}.
2310         @keyword upwind: value for upwind term M{upwind}
2311         @type upwind: any type that can be casted to L{Vector<escript.Vector>} object on L{Function<escript.Function>}.
2312         @keyword f: value for right hand side M{f}
2313         @type f: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
2314         @keyword alpha: value for right hand side M{S{alpha}}
2315         @type alpha: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
2316         @keyword g: value for right hand side M{g}
2317         @type g: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
2318         @keyword r: prescribed values M{r} for the solution in constraints.
2319         @type r: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2320                   depending of reduced order is used for the representation of the equation.
2321         @keyword q: mask for location of constraints
2322         @type q: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2323                   depending of reduced order is used for the representation of the equation.
2324         @raise IllegalCoefficient: if an unknown coefficient keyword is used.
2325         """
2326         super(AdvectionDiffusion, self).setValue(**coefficients)
2327    
2328       def getCoefficientOfGeneralPDE(self,name):
2329         """
2330         return the value of the coefficient name of the general PDE
2331    
2332         @param name: name of the coefficient requested.
2333         @type name: C{string}
2334         @return: the value of the coefficient  name
2335         @rtype: L{Data<escript.Data>}
2336         @raise IllegalCoefficient: if name is not one of coefficients
2337                      "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}.
2338         @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
2339         """
2340         if name == "A" :
2341             return self.getCoefficient("k")+util.outer(self.getCoefficient("upwind"),self.getCoefficient("upwind"))
2342         elif name == "B" :
2343             return escript.Data()
2344         elif name == "C" :
2345             return self.getCoefficient("v")
2346         elif name == "D" :
2347             return self.getCoefficient("omega")
2348         elif name == "X" :
2349             return escript.Data()
2350         elif name == "Y" :
2351             return self.getCoefficient("f")
2352         elif name == "d" :
2353             return self.getCoefficient("alpha")
2354         elif name == "y" :
2355             return self.getCoefficient("g")
2356         elif name == "d_contact" :
2357             return escript.Data()
2358         elif name == "y_contact" :
2359             return escript.Data()
2360         elif name == "r" :
2361             return self.getCoefficient("r")
2362         elif name == "q" :
2363             return self.getCoefficient("q")
2364         else:
2365            raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
2366    
2367    
2368  # $Log$  # $Log$
2369    # Revision 1.14  2005/09/22 01:54:57  jgs
2370    # Merge of development branch dev-02 back to main trunk on 2005-09-22
2371    #
2372    # Revision 1.13  2005/09/15 03:44:19  jgs
2373    # Merge of development branch dev-02 back to main trunk on 2005-09-15
2374    #
2375  # Revision 1.12  2005/09/01 03:31:28  jgs  # Revision 1.12  2005/09/01 03:31:28  jgs
2376  # Merge of development branch dev-02 back to main trunk on 2005-09-01  # Merge of development branch dev-02 back to main trunk on 2005-09-01
2377  #  #
# Line 2178  class AdvectivePDE(LinearPDE): Line 2381  class AdvectivePDE(LinearPDE):
2381  # Revision 1.10  2005/08/12 01:45:36  jgs  # Revision 1.10  2005/08/12 01:45:36  jgs
2382  # 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
2383  #  #
2384    # Revision 1.9.2.17  2005/09/21 07:03:33  matt
2385    # PDECoefficient and LinearPDE are now new style classes (introduced in Python
2386    # 2.2). Classes Poisson, Helmholtz, LameEquation and AdvectivePDE have been
2387    # modified to instead use portable/cooperative "super" calls to extend base
2388    # class methods.
2389    #
2390    # Revision 1.9.2.16  2005/09/16 01:54:37  matt
2391    # Removed redundant if-loop.
2392    #
2393    # Revision 1.9.2.15  2005/09/14 08:09:18  matt
2394    # Added "REDUCED" solution PDECoefficient descriptors for LinearPDEs.
2395    #
2396    # Revision 1.9.2.14  2005/09/07 06:26:16  gross
2397    # the solver from finley are put into the standalone package paso now
2398    #
2399  # Revision 1.9.2.13  2005/08/31 08:45:03  gross  # Revision 1.9.2.13  2005/08/31 08:45:03  gross
2400  # in the case of lumping no new system is allocated if the constraint is changed.  # in the case of lumping no new system is allocated if the constraint is changed.
2401  #  #

Legend:
Removed from v.149  
changed lines
  Added in v.425

  ViewVC Help
Powered by ViewVC 1.1.26