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

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

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

revision 328 by gross, Wed Dec 7 04:41:53 2005 UTC revision 425 by gross, Tue Jan 10 04:10:39 2006 UTC
# Line 17  the PDE solver library defined through t Line 17  the PDE solver library defined through t
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 449  class LinearPDE(object): Line 449  class LinearPDE(object):
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 498  class LinearPDE(object): Line 499  class LinearPDE(object):
499       self.__tolerance=1.e-8       self.__tolerance=1.e-8
500       self.__solver_method=self.DEFAULT       self.__solver_method=self.DEFAULT
501       self.__solver_package=self.DEFAULT       self.__solver_package=self.DEFAULT
502         self.__preconditioner=self.DEFAULT
503       self.__matrix_type=self.__domain.getSystemMatrixTypeId(self.DEFAULT,self.DEFAULT,False)       self.__matrix_type=self.__domain.getSystemMatrixTypeId(self.DEFAULT,self.DEFAULT,False)
504       self.__sym=False       self.__sym=False
505    
# Line 772  class LinearPDE(object): Line 774  class LinearPDE(object):
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 786  class LinearPDE(object): Line 786  class LinearPDE(object):
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()               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.")
# Line 815  class LinearPDE(object): Line 816  class LinearPDE(object):
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}, 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}.
825           @param preconditioner: sets a new solver method.
826           @type solver: one of L{DEFAULT}, L{JACOBI} L{ILU0}, L{ILUT},L{SSOR}
827         """         """
828         if solver==None: solve=self.DEFAULT         if solver==None: solve=self.DEFAULT
829         if not solver==self.getSolverMethod():         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 838  class LinearPDE(object): Line 843  class LinearPDE(object):
843    
844         m=self.getSolverMethod()         m=self.getSolverMethod()
845         p=self.getSolverPackage()         p=self.getSolverPackage()
846         if m==self.DEFAULT: method="DEFAULT"         method=""
847         elif m==self.DIRECT: method= "DIRECT"         if m[0]==self.DEFAULT: method="DEFAULT"
848         elif m==self.ITERATIVE: method= "ITERATIVE"         elif m[0]==self.DIRECT: method= "DIRECT"
849         elif m==self.CHOLEVSKY: method= "CHOLEVSKY"         elif m[0]==self.ITERATIVE: method= "ITERATIVE"
850         elif m==self.PCG: method= "PCG"         elif m[0]==self.CHOLEVSKY: method= "CHOLEVSKY"
851         elif m==self.CR: method= "CR"         elif m[0]==self.PCG: method= "PCG"
852         elif m==self.CGS: method= "CGS"         elif m[0]==self.CR: method= "CR"
853         elif m==self.BICGSTAB: method= "BICGSTAB"         elif m[0]==self.CGS: method= "CGS"
854         elif m==self.SSOR: method= "SSOR"         elif m[0]==self.BICGSTAB: method= "BICGSTAB"
855         elif m==self.GMRES: method= "GMRES"         elif m[0]==self.SSOR: method= "SSOR"
856         elif m==self.PRES20: method= "PRES20"         elif m[0]==self.GMRES: method= "GMRES"
857         elif m==self.LUMPING: method= "LUMPING"         elif m[0]==self.PRES20: method= "PRES20"
858         else : method="unknown"         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"         if p==self.DEFAULT: package="DEFAULT"
865         elif p==self.PASO: package= "PASO"         elif p==self.PASO: package= "PASO"
866         elif p==self.MKL: package= "MKL"         elif p==self.MKL: package= "MKL"
# Line 867  class LinearPDE(object): Line 877  class LinearPDE(object):
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):     def setSolverPackage(self,package=None):
883         """         """
# Line 898  class LinearPDE(object): Line 908  class LinearPDE(object):
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 1091  class LinearPDE(object): Line 1101  class LinearPDE(object):
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.getSolverPackage(),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 2242  class AdvectivePDE(LinearPDE): Line 2252  class AdvectivePDE(LinearPDE):
2252       elif name == "r" :       elif name == "r" :
2253           return self.getCoefficient("r")           return self.getCoefficient("r")
2254       elif name == "q" :       elif name == "q" :
2255             return self.getCoefficient("q")
2256         else:
2257            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")           return self.getCoefficient("q")
2364       else:       else:
2365          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
Legend:
Removed from v.328  
changed lines
  Added in v.425
  ViewVC Help
Powered by ViewVC 1.1.26