/[escript]/branches/symbolic_from_3470/escript/py_src/linearPDEs.py
ViewVC logotype

Diff of /branches/symbolic_from_3470/escript/py_src/linearPDEs.py

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

revision 2473 by gross, Thu Jun 11 08:32:32 2009 UTC revision 2474 by gross, Tue Jun 16 06:32:15 2009 UTC
# Line 172  class SolverOptions(object): Line 172  class SolverOptions(object):
172          out+="\nAbsolute tolerance = %e"%self.getAbsoluteTolerance()          out+="\nAbsolute tolerance = %e"%self.getAbsoluteTolerance()
173          out+="\nSymmetric problem = %s"%self.isSymmetric()          out+="\nSymmetric problem = %s"%self.isSymmetric()
174          out+="\nMaximum number of iteration steps = %s"%self.getIterMax()          out+="\nMaximum number of iteration steps = %s"%self.getIterMax()
175            out+="\nInner tolerance = %e"%self.getInnerTolerance()
176            out+="\nAdapt innner tolerance = %s"%self.adaptInnerTolerance()
177        
178          if self.getPackage() == self.PASO:          if self.getPackage() == self.PASO:
179              out+="\nSolver method = %s"%self.getName(self.getSolverMethod())              out+="\nSolver method = %s"%self.getName(self.getSolverMethod())
180              if self.getSolverMethod() == self.GMRES:              if self.getSolverMethod() == self.GMRES:
# Line 331  class SolverOptions(object): Line 334  class SolverOptions(object):
334          Sets the key of the coarsening method to be applied in AMG.          Sets the key of the coarsening method to be applied in AMG.
335    
336          @param method: selects the coarsening method .          @param method: selects the coarsening method .
337          @type method: in L{SolverOptions.DEFAULT}, L{SolverOptions.YAIR_SHAPIRA_COARSENING},          @type method: in {SolverOptions.DEFAULT}, L{SolverOptions.YAIR_SHAPIRA_COARSENING},
338          L{SolverOptions.RUGE_STUEBEN_COARSENING}, L{SolverOptions.AGGREGATION_COARSENING}          L{SolverOptions.RUGE_STUEBEN_COARSENING}, L{SolverOptions.AGGREGATION_COARSENING}
339          """          """
340        if method==None: method=0
341          if not method in [self.DEFAULT, self.YAIR_SHAPIRA_COARSENING, self.RUGE_STUEBEN_COARSENING, self.AGGREGATION_COARSENING]:          if not method in [self.DEFAULT, self.YAIR_SHAPIRA_COARSENING, self.RUGE_STUEBEN_COARSENING, self.AGGREGATION_COARSENING]:
342               raise ValueError,"unknown coarsening method %s"%method               raise ValueError,"unknown coarsening method %s"%method
343          self.__coarsening=method          self.__coarsening=method
# Line 356  class SolverOptions(object): Line 360  class SolverOptions(object):
360          @note: Not all packages support all preconditioner. It can be assumed that a package makes a reasonable choice if it encounters          @note: Not all packages support all preconditioner. It can be assumed that a package makes a reasonable choice if it encounters
361          an unknown preconditioner.          an unknown preconditioner.
362          """          """
363        if preconditioner==None: preconditioner=10
364          if not preconditioner in [ SolverOptions.SSOR, SolverOptions.ILU0, SolverOptions.ILUT, SolverOptions.JACOBI,          if not preconditioner in [ SolverOptions.SSOR, SolverOptions.ILU0, SolverOptions.ILUT, SolverOptions.JACOBI,
365                                      SolverOptions.AMG, SolverOptions.REC_ILU, SolverOptions.GAUSS_SEIDEL, SolverOptions.RILU,                                      SolverOptions.AMG, SolverOptions.REC_ILU, SolverOptions.GAUSS_SEIDEL, SolverOptions.RILU,
366                                      SolverOptions.NO_PRECONDITIONER] :                                      SolverOptions.NO_PRECONDITIONER] :
# Line 385  class SolverOptions(object): Line 390  class SolverOptions(object):
390          @note: Not all packages support all solvers. It can be assumed that a package makes a reasonable choice if it encounters          @note: Not all packages support all solvers. It can be assumed that a package makes a reasonable choice if it encounters
391          an unknown solver method.          an unknown solver method.
392          """          """
393        if method==None: method=0
394          if not method in [ SolverOptions.DEFAULT, SolverOptions.DIRECT, SolverOptions.CHOLEVSKY, SolverOptions.PCG,          if not method in [ SolverOptions.DEFAULT, SolverOptions.DIRECT, SolverOptions.CHOLEVSKY, SolverOptions.PCG,
395                             SolverOptions.CR, SolverOptions.CGS, SolverOptions.BICGSTAB, SolverOptions.SSOR,                             SolverOptions.CR, SolverOptions.CGS, SolverOptions.BICGSTAB, SolverOptions.SSOR,
396                             SolverOptions.GMRES, SolverOptions.PRES20, SolverOptions.LUMPING, SolverOptions.ITERATIVE, SolverOptions.AMG,                             SolverOptions.GMRES, SolverOptions.PRES20, SolverOptions.LUMPING, SolverOptions.ITERATIVE, SolverOptions.AMG,
# Line 411  class SolverOptions(object): Line 417  class SolverOptions(object):
417          @type package: in L{SolverOptions.DEFAULT}, L{SolverOptions.PASO}, L{SolverOptions.SUPER_LU}, L{SolverOptions.PASTIX}, L{SolverOptions.MKL}, L{SolverOptions.UMFPACK}, L{SolverOptions.TRILINOS}          @type package: in L{SolverOptions.DEFAULT}, L{SolverOptions.PASO}, L{SolverOptions.SUPER_LU}, L{SolverOptions.PASTIX}, L{SolverOptions.MKL}, L{SolverOptions.UMFPACK}, L{SolverOptions.TRILINOS}
418          @note: Not all packages are support on all implementation. An exception may be thrown on some platforms if a particular is requested.          @note: Not all packages are support on all implementation. An exception may be thrown on some platforms if a particular is requested.
419          """          """
420        if package==None: package=0
421          if not package in [SolverOptions.DEFAULT, SolverOptions.PASO, SolverOptions.SUPER_LU, SolverOptions.PASTIX, SolverOptions.MKL, SolverOptions.UMFPACK, SolverOptions.TRILINOS]:          if not package in [SolverOptions.DEFAULT, SolverOptions.PASO, SolverOptions.SUPER_LU, SolverOptions.PASTIX, SolverOptions.MKL, SolverOptions.UMFPACK, SolverOptions.TRILINOS]:
422               raise ValueError,"unknown solver package %s"%package               raise ValueError,"unknown solver package %s"%package
423          self.__package=package          self.__package=package
# Line 456  class SolverOptions(object): Line 463  class SolverOptions(object):
463              if restart<1:              if restart<1:
464                  raise ValueError,"restart must be positive."                  raise ValueError,"restart must be positive."
465              self.__restart=restart              self.__restart=restart
466            
467      def getRestart(self):      def getRestart(self):
468          """          """
469          Returns the number of iterations steps after which GMRES is performing a restart.          Returns the number of iterations steps after which GMRES is performing a restart.
# Line 467  class SolverOptions(object): Line 475  class SolverOptions(object):
475              return None              return None
476          else:          else:
477              return self.__restart              return self.__restart
478        def _getRestartForC(self):
479            r=self.getRestart()
480            if r == None:
481                return -1
482                else:
483                return r
484      def setTruncation(self,truncation=20):      def setTruncation(self,truncation=20):
485          """          """
486          Sets the number of residuals in GMRES to be stored for orthogonalization.  The more residuals are stored          Sets the number of residuals in GMRES to be stored for orthogonalization.  The more residuals are stored
# Line 775  class SolverOptions(object): Line 789  class SolverOptions(object):
789          Switches the verbosity of the solver off.          Switches the verbosity of the solver off.
790          """          """
791          self.__verbose=False          self.__verbose=False
792      def setVerbosity(self,verbose=True):      def setVerbosity(self,verbose=False):
793          """          """
794          Sets the verbosity flag for the solver to C{flag}.          Sets the verbosity flag for the solver to C{flag}.
795    
# Line 786  class SolverOptions(object): Line 800  class SolverOptions(object):
800              self.setVerbosityOn()              self.setVerbosityOn()
801          else:          else:
802              self.setVerbosityOff()              self.setVerbosityOff()
803            
         self.__adapt_inner_tolerance=True  
804      def adaptInnerTolerance(self):      def adaptInnerTolerance(self):
805          """          """
806          Returns C{True} if the tolerance of the inner solver is selected automatically.          Returns C{True} if the tolerance of the inner solver is selected automatically.
# Line 1208  class LinearProblem(object): Line 1221  class LinearProblem(object):
1221     where M{L} is an operator and M{f} is the right hand side. This operator     where M{L} is an operator and M{f} is the right hand side. This operator
1222     problem will be solved to get the unknown M{u}.     problem will be solved to get the unknown M{u}.
1223    
1224     @cvar DEFAULT: The default method used to solve the system of linear     """
                   equations  
    @cvar DIRECT: The direct solver based on LDU factorization  
    @cvar CHOLEVSKY: The direct solver based on LDLt factorization (can only be  
                     applied for symmetric PDEs)  
    @cvar PCG: The preconditioned conjugate gradient method (can only be applied  
               for symmetric PDEs)  
    @cvar CR: The conjugate residual method  
    @cvar CGS: The conjugate gradient square method  
    @cvar BICGSTAB: The stabilized BiConjugate Gradient method  
    @cvar TFQMR: Transport Free Quasi Minimal Residual method  
    @cvar MINRES: Minimum residual method  
    @cvar SSOR: The symmetric overrelaxation method  
    @cvar ILU0: The incomplete LU factorization preconditioner with no fill-in  
    @cvar ILUT: The incomplete LU factorization preconditioner with fill-in  
    @cvar JACOBI: The Jacobi preconditioner  
    @cvar GMRES: The Gram-Schmidt minimum residual method  
    @cvar PRES20: Special GMRES with restart after 20 steps and truncation after  
                  5 residuals  
    @cvar LUMPING: Matrix lumping  
    @cvar NO_REORDERING: No matrix reordering allowed  
    @cvar MINIMUM_FILL_IN: Reorder matrix to reduce fill-in during factorization  
    @cvar NESTED_DISSECTION: Reorder matrix to improve load balancing during  
                             factorization  
    @cvar PASO: PASO solver package  
    @cvar SCSL: SGI SCSL solver library  
    @cvar MKL: Intel's MKL solver library  
    @cvar UMFPACK: The UMFPACK library  
    @cvar TRILINOS: The TRILINOS parallel solver class library from Sandia Natl  
                    Labs  
    @cvar ITERATIVE: The default iterative solver  
    @cvar AMG: Algebraic Multi Grid  
    @cvar RILU: Recursive ILU  
    @cvar GS: Gauss-Seidel solver  
   
    """  
    DEFAULT= 0  
    DIRECT= 1  
    CHOLEVSKY= 2  
    PCG= 3  
    CR= 4  
    CGS= 5  
    BICGSTAB= 6  
    SSOR= 7  
    ILU0= 8  
    ILUT= 9  
    JACOBI= 10  
    GMRES= 11  
    PRES20= 12  
    LUMPING= 13  
    NO_REORDERING= 17  
    MINIMUM_FILL_IN= 18  
    NESTED_DISSECTION= 19  
    SCSL= 14  
    MKL= 15  
    UMFPACK= 16  
    ITERATIVE= 20  
    PASO= 21  
    AMG= 22  
    RILU = 23  
    TRILINOS = 24  
    NONLINEAR_GMRES = 25  
    TFQMR = 26  
    MINRES = 27  
    GS=28  
   
    SMALL_TOLERANCE=1.e-13  
    PACKAGE_KEY="package"  
    METHOD_KEY="method"  
    SYMMETRY_KEY="symmetric"  
    TOLERANCE_KEY="tolerance"  
    PRECONDITIONER_KEY="preconditioner"  
   
   
1225     def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):     def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):
1226       """       """
1227       Initializes a linear problem.       Initializes a linear problem.
# Line 1305  class LinearProblem(object): Line 1245  class LinearProblem(object):
1245       self.__altered_coefficients=False       self.__altered_coefficients=False
1246       self.__reduce_equation_order=False       self.__reduce_equation_order=False
1247       self.__reduce_solution_order=False       self.__reduce_solution_order=False
1248       self.__tolerance=1.e-8       self.__sym=False
1249       self.__solver_method=self.DEFAULT       self.setSolverOptions()
      self.__solver_package=self.DEFAULT  
      self.__preconditioner=self.DEFAULT  
1250       self.__is_RHS_valid=False       self.__is_RHS_valid=False
1251       self.__is_operator_valid=False       self.__is_operator_valid=False
      self.__sym=False  
1252       self.__COEFFICIENTS={}       self.__COEFFICIENTS={}
1253         self.__solution_rtol=1.e99
1254         self.__solution_atol=1.e99
1255         self.setSymmetryOff()
1256       # initialize things:       # initialize things:
1257       self.resetAllCoefficients()       self.resetAllCoefficients()
1258       self.__system_type=None       self.initializeSystem()
      self.updateSystemType()  
1259     # ==========================================================================     # ==========================================================================
1260     #    general stuff:     #    general stuff:
1261     # ==========================================================================     # ==========================================================================
# Line 1481  class LinearProblem(object): Line 1420  class LinearProblem(object):
1420     # ==========================================================================     # ==========================================================================
1421     #   solver settings:     #   solver settings:
1422     # ==========================================================================     # ==========================================================================
1423     def setSolverMethod(self,solver=None,preconditioner=None):     def setSolverOptions(self,options=None):
        """  
        Sets a new solver method and/or preconditioner.  
   
        @param solver: new solver method to use  
        @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{TFQMR}, L{MINRES}, L{PRES20}, L{LUMPING}, L{AMG}  
        @param preconditioner: new preconditioner to use  
        @type preconditioner: one of L{DEFAULT}, L{JACOBI} L{ILU0}, L{ILUT},  
                              L{SSOR}, L{RILU}, L{GS}  
   
        @note: the solver method chosen may not be available by the selected  
               package.  
        """  
        if solver==None: solver=self.__solver_method  
        if preconditioner==None: preconditioner=self.__preconditioner  
        if solver==None: solver=self.DEFAULT  
        if preconditioner==None: preconditioner=self.DEFAULT  
        if not (solver,preconditioner)==self.getSolverMethod():  
            self.__solver_method=solver  
            self.__preconditioner=preconditioner  
            self.updateSystemType()  
            self.trace("New solver is %s"%self.getSolverMethodName())  
   
    def getSolverMethodName(self):  
        """  
        Returns the name of the solver currently used.  
   
        @return: the name of the solver currently used  
        @rtype: C{string}  
        """  
   
        m=self.getSolverMethod()  
        p=self.getSolverPackage()  
        method=""  
        if m[0]==self.DEFAULT: method="DEFAULT"  
        elif m[0]==self.DIRECT: method= "DIRECT"  
        elif m[0]==self.ITERATIVE: method= "ITERATIVE"  
        elif m[0]==self.CHOLEVSKY: method= "CHOLEVSKY"  
        elif m[0]==self.PCG: method= "PCG"  
        elif m[0]==self.TFQMR: method= "TFQMR"  
        elif m[0]==self.MINRES: method= "MINRES"  
        elif m[0]==self.CR: method= "CR"  
        elif m[0]==self.CGS: method= "CGS"  
        elif m[0]==self.BICGSTAB: method= "BICGSTAB"  
        elif m[0]==self.SSOR: method= "SSOR"  
        elif m[0]==self.GMRES: method= "GMRES"  
        elif m[0]==self.PRES20: method= "PRES20"  
        elif m[0]==self.LUMPING: method= "LUMPING"  
        elif m[0]==self.AMG: method= "AMG"  
        if m[1]==self.DEFAULT: method+="+DEFAULT"  
        elif m[1]==self.JACOBI: method+= "+JACOBI"  
        elif m[1]==self.ILU0: method+= "+ILU0"  
        elif m[1]==self.ILUT: method+= "+ILUT"  
        elif m[1]==self.SSOR: method+= "+SSOR"  
        elif m[1]==self.AMG: method+= "+AMG"  
        elif m[1]==self.RILU: method+= "+RILU"  
        elif m[1]==self.GS: method+= "+GS"  
        if p==self.DEFAULT: package="DEFAULT"  
        elif p==self.PASO: package= "PASO"  
        elif p==self.MKL: package= "MKL"  
        elif p==self.SCSL: package= "SCSL"  
        elif p==self.UMFPACK: package= "UMFPACK"  
        elif p==self.TRILINOS: package= "TRILINOS"  
        else : method="unknown"  
        return "%s solver of %s package"%(method,package)  
   
    def getSolverMethod(self):  
        """  
        Returns the solver method and preconditioner used.  
   
        @return: the solver method currently used.  
        @rtype: C{tuple} of C{int}  
        """  
        return self.__solver_method,self.__preconditioner  
   
    def setSolverPackage(self,package=None):  
        """  
        Sets a new solver package.  
   
        @param package: the new solver package  
        @type package: one of L{DEFAULT}, L{PASO} L{SCSL}, L{MKL}, L{UMFPACK},  
                       L{TRILINOS}  
        """  
        if package==None: package=self.DEFAULT  
        if not package==self.getSolverPackage():  
            self.__solver_package=package  
            self.updateSystemType()  
            self.trace("New solver is %s"%self.getSolverMethodName())  
   
    def getSolverPackage(self):  
        """  
        Returns the package of the solver.  
   
        @return: the solver package currently being used  
        @rtype: C{int}  
1424         """         """
1425         return self.__solver_package         Sets the solver options.
1426    
1427           @param options: the new solver options. If equal C{None}, the solver options are set to the default.
1428           @type options: L{SolverOptions} or C{None}
1429           @note: The symmetry flag of options is overwritten by the symmetry flag of the L{LinearProblem}.
1430           """
1431           if options==None:
1432              self.__solver_options=SolverOptions()
1433           elif isinstance(options, SolverOptions):
1434              self.__solver_options=options
1435           else:
1436              raise ValueError,"options must be a SolverOptions object."
1437           self.__solver_options.setSymmetry(self.__sym)
1438        
1439       def getSolverOptions(self):
1440           """
1441           Returns the solver options
1442      
1443           @rtype: L{SolverOptions}
1444           """
1445           self.__solver_options.setSymmetry(self.__sym)
1446           return self.__solver_options
1447          
1448     def isUsingLumping(self):     def isUsingLumping(self):
1449        """        """
1450        Checks if matrix lumping is the current solver method.        Checks if matrix lumping is the current solver method.
# Line 1588  class LinearProblem(object): Line 1452  class LinearProblem(object):
1452        @return: True if the current solver method is lumping        @return: True if the current solver method is lumping
1453        @rtype: C{bool}        @rtype: C{bool}
1454        """        """
1455        return self.getSolverMethod()[0]==self.LUMPING        return self.getSolverOptions().getSolverMethod()==self.getSolverOptions().LUMPING
   
    def setTolerance(self,tol=1.e-8):  
        """  
        Resets the tolerance for the solver method to C{tol}.  
   
        @param tol: new tolerance for the solver. If C{tol} is lower than the  
                    current tolerence the system will be resolved.  
        @type tol: positive C{float}  
        @raise ValueError: if tolerance is not positive  
        """  
        if not tol>0:  
            raise ValueError,"Tolerance has to be positive"  
        if tol<self.getTolerance(): self.invalidateSolution()  
        self.trace("New tolerance %e"%tol)  
        self.__tolerance=tol  
        return  
   
    def getTolerance(self):  
        """  
        Returns the tolerance currently set for the solution.  
   
        @return: tolerance currently used  
        @rtype: C{float}  
        """  
        return self.__tolerance  
   
1456     # ==========================================================================     # ==========================================================================
1457     #    symmetry  flag:     #    symmetry  flag:
1458     # ==========================================================================     # ==========================================================================
# Line 1624  class LinearProblem(object): Line 1462  class LinearProblem(object):
1462    
1463        @return: True if a symmetric PDE is indicated, False otherwise        @return: True if a symmetric PDE is indicated, False otherwise
1464        @rtype: C{bool}        @rtype: C{bool}
1465          @note: the method is equivalent to use getSolverOptions().isSymmetric()
1466        """        """
1467        return self.__sym        self.getSolverOptions().isSymmetric()
1468    
1469     def setSymmetryOn(self):     def setSymmetryOn(self):
1470        """        """
1471        Sets the symmetry flag.        Sets the symmetry flag.
1472          @note: The method overwrites the symmetry flag set by the solver options
1473        """        """
1474        if not self.isSymmetric():        self.__sym=True
1475           self.trace("PDE is set to be symmetric")        self.getSolverOptions().setSymmetryOn()
          self.__sym=True  
          self.updateSystemType()  
1476    
1477     def setSymmetryOff(self):     def setSymmetryOff(self):
1478        """        """
1479        Clears the symmetry flag.        Clears the symmetry flag.
1480          @note: The method overwrites the symmetry flag set by the solver options
1481        """        """
1482        if self.isSymmetric():        self.__sym=False
1483           self.trace("PDE is set to be nonsymmetric")        self.getSolverOptions().setSymmetryOff()
          self.__sym=False  
          self.updateSystemType()  
1484    
1485     def setSymmetryTo(self,flag=False):     def setSymmetry(self,flag=False):
1486        """        """
1487        Sets the symmetry flag to C{flag}.        Sets the symmetry flag to C{flag}.
1488    
1489        @param flag: If True, the symmetry flag is set otherwise reset.        @param flag: If True, the symmetry flag is set otherwise reset.
1490        @type flag: C{bool}        @type flag: C{bool}
1491          @note: The method overwrites the symmetry flag set by the solver options
1492        """        """
1493        if flag:        self.getSolverOptions().setSymmetry(flag)
          self.setSymmetryOn()  
       else:  
          self.setSymmetryOff()  
   
1494     # ==========================================================================     # ==========================================================================
1495     # function space handling for the equation as well as the solution     # function space handling for the equation as well as the solution
1496     # ==========================================================================     # ==========================================================================
# Line 1783  class LinearProblem(object): Line 1617  class LinearProblem(object):
1617          self.setReducedOrderForEquationOn()          self.setReducedOrderForEquationOn()
1618       else:       else:
1619          self.setReducedOrderForEquationOff()          self.setReducedOrderForEquationOff()
1620       def getOperatorType(self):
    def updateSystemType(self):  
      """  
      Reassesses the system type and, if a new matrix is needed, resets the  
      system.  
      """  
      new_system_type=self.getRequiredSystemType()  
      if not new_system_type==self.__system_type:  
          self.trace("Matrix type is now %d."%new_system_type)  
          self.__system_type=new_system_type  
          self.initializeSystem()  
   
    def getSystemType(self):  
1621        """        """
1622        Returns the current system type.        Returns the current system type.
1623        """        """
1624        return self.__system_type        return self.__operator_type
1625    
1626     def checkSymmetricTensor(self,name,verbose=True):     def checkSymmetricTensor(self,name,verbose=True):
1627        """        """
# Line 1813  class LinearProblem(object): Line 1635  class LinearProblem(object):
1635        @return: True if coefficient C{name} is symmetric        @return: True if coefficient C{name} is symmetric
1636        @rtype: C{bool}        @rtype: C{bool}
1637        """        """
1638          SMALL_TOLERANCE=util.EPSILON*10.
1639        A=self.getCoefficient(name)        A=self.getCoefficient(name)
1640        verbose=verbose or self.__debug        verbose=verbose or self.__debug
1641        out=True        out=True
1642        if not A.isEmpty():        if not A.isEmpty():
1643           tol=util.Lsup(A)*self.SMALL_TOLERANCE           tol=util.Lsup(A)*SMALL_TOLERANCE
1644           s=A.getShape()           s=A.getShape()
1645           if A.getRank() == 4:           if A.getRank() == 4:
1646              if s[0]==s[2] and s[1] == s[3]:              if s[0]==s[2] and s[1] == s[3]:
# Line 1862  class LinearProblem(object): Line 1685  class LinearProblem(object):
1685                 symmetric.                 symmetric.
1686        @rtype: C{bool}        @rtype: C{bool}
1687        """        """
1688          SMALL_TOLERANCE=util.EPSILON*10.
1689        B=self.getCoefficient(name0)        B=self.getCoefficient(name0)
1690        C=self.getCoefficient(name1)        C=self.getCoefficient(name1)
1691        verbose=verbose or self.__debug        verbose=verbose or self.__debug
# Line 1875  class LinearProblem(object): Line 1699  class LinearProblem(object):
1699        elif not B.isEmpty() and not C.isEmpty():        elif not B.isEmpty() and not C.isEmpty():
1700           sB=B.getShape()           sB=B.getShape()
1701           sC=C.getShape()           sC=C.getShape()
1702           tol=(util.Lsup(B)+util.Lsup(C))*self.SMALL_TOLERANCE/2.           tol=(util.Lsup(B)+util.Lsup(C))*SMALL_TOLERANCE/2.
1703           if len(sB) != len(sC):           if len(sB) != len(sC):
1704               if verbose: print "non-symmetric problem because ranks of %s (=%s) and %s (=%s) are different."%(name0,len(sB),name1,len(sC))               if verbose: print "non-symmetric problem because ranks of %s (=%s) and %s (=%s) are different."%(name0,len(sB),name1,len(sC))
1705               out=False               out=False
# Line 2020  class LinearProblem(object): Line 1844  class LinearProblem(object):
1844         """         """
1845         Returns True if the solution is still valid.         Returns True if the solution is still valid.
1846         """         """
1847           if self.__solution_rtol>self.getSolverOptions().getTolerance() or \
1848              self.__solution_atol>self.getSolverOptions().getAbsoluteTolerance():
1849             self.invalidateSolution()  
1850         return self.__is_solution_valid         return self.__is_solution_valid
1851    
1852     def validOperator(self):     def validOperator(self):
# Line 2040  class LinearProblem(object): Line 1867  class LinearProblem(object):
1867         """         """
1868         Returns True if the operator is still valid.         Returns True if the operator is still valid.
1869         """         """
1870           if self.getRequiredOperatorType()==self.getOperatorType(): self.invalidateOperator()
1871         return self.__is_operator_valid         return self.__is_operator_valid
1872    
1873     def validRightHandSide(self):     def validRightHandSide(self):
# Line 2066  class LinearProblem(object): Line 1894  class LinearProblem(object):
1894         """         """
1895         Announces that everything has to be rebuilt.         Announces that everything has to be rebuilt.
1896         """         """
        if self.isRightHandSideValid(): self.trace("System has to be rebuilt.")  
1897         self.invalidateSolution()         self.invalidateSolution()
1898         self.invalidateOperator()         self.invalidateOperator()
1899         self.invalidateRightHandSide()         self.invalidateRightHandSide()
# Line 2082  class LinearProblem(object): Line 1909  class LinearProblem(object):
1909         Resets the system clearing the operator, right hand side and solution.         Resets the system clearing the operator, right hand side and solution.
1910         """         """
1911         self.trace("New System has been created.")         self.trace("New System has been created.")
1912           self.__operator_type=None
1913         self.__operator=escript.Operator()         self.__operator=escript.Operator()
1914         self.__righthandside=escript.Data()         self.__righthandside=escript.Data()
1915         self.__solution=escript.Data()         self.__solution=escript.Data()
# Line 2136  class LinearProblem(object): Line 1964  class LinearProblem(object):
1964    
1965     def setSolution(self,u):     def setSolution(self,u):
1966         """         """
1967         Sets the solution assuming that makes the system valid.         Sets the solution assuming that makes the system valid with the tolrance
1968           defined by the solver options
1969         """         """
1970           self.__solution_rtol=self.getSolverOptions().getTolerance()
1971           self.__solution_atol=self.getSolverOptions().getAbsoluteTolerance()
1972         self.__solution=u         self.__solution=u
1973         self.validSolution()         self.validSolution()
1974    
# Line 2170  class LinearProblem(object): Line 2001  class LinearProblem(object):
2001         Makes sure that the operator is instantiated and returns it initialized         Makes sure that the operator is instantiated and returns it initialized
2002         with zeros.         with zeros.
2003         """         """
2004         if self.__operator.isEmpty():         if self.getOperatorType() == None:
2005             if self.isUsingLumping():             if self.isUsingLumping():
2006                 self.__operator=self.createSolution()                 self.__operator=self.createSolution()
2007             else:             else:
2008                 self.__operator=self.createOperator()                 self.__operator=self.createOperator()
2009           self.__operator_type=self.getRequiredOperatorType()
2010         else:         else:
2011             if self.isUsingLumping():             if self.isUsingLumping():
2012                 self.__operator.setToZero()                 self.__operator.setToZero()
# Line 2186  class LinearProblem(object): Line 2018  class LinearProblem(object):
2018         """         """
2019         Returns the operator in its current state.         Returns the operator in its current state.
2020         """         """
        if self.__operator.isEmpty(): self.resetOperator()  
2021         return self.__operator         return self.__operator
2022    
2023     def setValue(self,**coefficients):     def setValue(self,**coefficients):
# Line 2248  class LinearProblem(object): Line 2079  class LinearProblem(object):
2079     # methods that are typically overwritten when implementing a particular     # methods that are typically overwritten when implementing a particular
2080     # linear problem     # linear problem
2081     # ==========================================================================     # ==========================================================================
2082     def getRequiredSystemType(self):     def getRequiredOperatorType(self):
2083        """        """
2084        Returns the system type which needs to be used by the current set up.        Returns the system type which needs to be used by the current set up.
2085    
2086        @note: Typically this method is overwritten when implementing a        @note: Typically this method is overwritten when implementing a
2087               particular linear problem.               particular linear problem.
2088        """        """
2089        return 0        return None
2090    
2091     def createOperator(self):     def createOperator(self):
2092         """         """
# Line 2304  class LinearProblem(object): Line 2135  class LinearProblem(object):
2135                linear problem.                linear problem.
2136         """         """
2137         return (self.getCurrentOperator(), self.getCurrentRightHandSide())         return (self.getCurrentOperator(), self.getCurrentRightHandSide())
 #=====================  
2138    
2139  class LinearPDE(LinearProblem):  class LinearPDE(LinearProblem):
2140     """     """
# Line 2489  class LinearPDE(LinearProblem): Line 2319  class LinearPDE(LinearProblem):
2319       """       """
2320       return "<LinearPDE %d>"%id(self)       return "<LinearPDE %d>"%id(self)
2321    
2322     def getRequiredSystemType(self):     def getRequiredOperatorType(self):
2323        """        """
2324        Returns the system type which needs to be used by the current set up.        Returns the system type which needs to be used by the current set up.
2325        """        """
2326        return self.getDomain().getSystemMatrixTypeId(self.getSolverMethod()[0],self.getSolverMethod()[1],self.getSolverPackage(),self.isSymmetric())        solver_options=self.getSolverOptions()
2327          return self.getDomain().getSystemMatrixTypeId(solver_options.getSolverMethod(), solver_options.getPreconditioner(),solver_options.getPackage(), solver_options.isSymmetric())
2328    
2329     def checkSymmetry(self,verbose=True):     def checkSymmetry(self,verbose=True):
2330        """        """
# Line 2524  class LinearPDE(LinearProblem): Line 2355  class LinearPDE(LinearProblem):
2355         """         """
2356         Returns an instance of a new operator.         Returns an instance of a new operator.
2357         """         """
2358         self.trace("New operator is allocated.")         optype=self.getRequiredOperatorType()
2359           self.trace("New operator of type %s is allocated."%optype)
2360         return self.getDomain().newOperator( \         return self.getDomain().newOperator( \
2361                             self.getNumEquations(), \                             self.getNumEquations(), \
2362                             self.getFunctionSpaceForEquation(), \                             self.getFunctionSpaceForEquation(), \
2363                             self.getNumSolutions(), \                             self.getNumSolutions(), \
2364                             self.getFunctionSpaceForSolution(), \                             self.getFunctionSpaceForSolution(), \
2365                             self.getSystemType())                             optype)
2366    
2367     def getSolution(self,**options):     def getSolution(self):
2368         """         """
2369         Returns the solution of the PDE.         Returns the solution of the PDE.
2370    
2371         @return: the solution         @return: the solution
2372         @rtype: L{Data<escript.Data>}         @rtype: L{Data<escript.Data>}
        @param options: solver options  
        @keyword verbose: True to get some information during PDE solution  
        @type verbose: C{bool}  
        @keyword reordering: reordering scheme to be used during elimination.  
                             Allowed values are L{NO_REORDERING},  
                             L{MINIMUM_FILL_IN} and L{NESTED_DISSECTION}  
        @keyword iter_max: maximum number of iteration steps allowed  
        @keyword drop_tolerance: threshold for dropping in L{ILUT}  
        @keyword drop_storage: maximum of allowed memory in L{ILUT}  
        @keyword truncation: maximum number of residuals in L{GMRES}  
        @keyword restart: restart cycle length in L{GMRES}  
2373         """         """
2374           option_class=self.getSolverOptions()
2375         if not self.isSolutionValid():         if not self.isSolutionValid():
2376            mat,f=self.getSystem()            mat,f=self.getSystem()
2377            if self.isUsingLumping():            if self.isUsingLumping():
2378               self.setSolution(f*1/mat)               self.setSolution(f*1/mat)
2379            else:            else:
              options[self.TOLERANCE_KEY]=self.getTolerance()  
              options[self.METHOD_KEY]=self.getSolverMethod()[0]  
              options[self.PRECONDITIONER_KEY]=self.getSolverMethod()[1]  
              options[self.PACKAGE_KEY]=self.getSolverPackage()  
              options[self.SYMMETRY_KEY]=self.isSymmetric()  
2380               self.trace("PDE is resolved.")               self.trace("PDE is resolved.")
2381               self.trace("solver options: %s"%str(options))               self.trace("solver options: %s"%str(option_class))
2382               self.setSolution(mat.solve(f,options))               self.setSolution(mat.solve(f,option_class))
2383         return self.getCurrentSolution()         return self.getCurrentSolution()
2384    
2385     def getSystem(self):     def getSystem(self):
# Line 3535  class TransportPDE(LinearProblem): Line 3352  class TransportPDE(LinearProblem):
3352           theta=1.           theta=1.
3353         else:         else:
3354           theta=0.5           theta=0.5
3355           optype=self.getRequiredOperatorType()
3356         self.trace("New Transport problem is allocated.")         self.trace("New Transport problem pf type %s is allocated."%optype)
3357         return self.getDomain().newTransportProblem( \         return self.getDomain().newTransportProblem( \
3358                                 theta,                                 theta,
3359                                 self.getNumEquations(), \                                 self.getNumEquations(), \
3360                                 self.getFunctionSpaceForSolution(), \                                 self.getFunctionSpaceForSolution(), \
3361                                 self.getSystemType())                                 optype)
3362    
3363     def setInitialSolution(self,u):     def setInitialSolution(self,u):
3364         """         """
# Line 3561  class TransportPDE(LinearProblem): Line 3378  class TransportPDE(LinearProblem):
3378                raise ValueError,"Illegal shape %s of initial solution."%(u2.getShape(),)                raise ValueError,"Illegal shape %s of initial solution."%(u2.getShape(),)
3379         self.getOperator().setInitialValue(u2)         self.getOperator().setInitialValue(u2)
3380    
3381     def getRequiredSystemType(self):     def getRequiredOperatorType(self):
3382        """        """
3383        Returns the system type which needs to be used by the current set up.        Returns the system type which needs to be used by the current set up.
3384    
3385        @return: a code to indicate the type of transport problem scheme used        @return: a code to indicate the type of transport problem scheme used
3386        @rtype: C{float}        @rtype: C{float}
3387        """        """
3388        return self.getDomain().getTransportTypeId(self.getSolverMethod()[0],self.getSolverMethod()[1],self.getSolverPackage(),self.isSymmetric())        solver_options=self.getSolverOptions()
3389          return self.getDomain().getTransportTypeId(solver_options.getSolverMethod(), solver_options.getPreconditioner(),solver_options.getPackage(), solver_options.isSymmetric())
3390    
3391     def getUnlimitedTimeStepSize(self):     def getUnlimitedTimeStepSize(self):
3392        """        """
# Line 3613  class TransportPDE(LinearProblem): Line 3431  class TransportPDE(LinearProblem):
3431         """         """
3432         return self.__constraint_factor         return self.__constraint_factor
3433     #====================================================================     #====================================================================
3434     def getSolution(self,dt,**options):     def getSolution(self,dt):
3435         """         """
3436         Returns the solution of the problem.         Returns the solution of the problem.
3437    
3438         @return: the solution         @return: the solution
3439         @rtype: L{Data<escript.Data>}         @rtype: L{Data<escript.Data>}
   
3440         """         """
3441           option_class=self.getSolverOptions()
3442         if dt<=0:         if dt<=0:
3443             raise ValueError,"step size needs to be positive."             raise ValueError,"step size needs to be positive."
3444         options[self.TOLERANCE_KEY]=self.getTolerance()         self.setSolution(self.getOperator().solve(self.getRightHandSide(),dt,option_class))
        self.setSolution(self.getOperator().solve(self.getRightHandSide(),dt,options))  
3445         self.validSolution()         self.validSolution()
3446         return self.getCurrentSolution()         return self.getCurrentSolution()
3447    
Legend:
Removed from v.2473  
changed lines
  Added in v.2474
  ViewVC Help
Powered by ViewVC 1.1.26