/[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 2470 by gross, Thu Jun 11 08:32:32 2009 UTC revision 2549 by jfenwick, Mon Jul 20 06:43:47 2009 UTC
# Line 1  Line 1 
1    
2  ########################################################  ########################################################
3  #  #
4  # Copyright (c) 2003-2008 by University of Queensland  # Copyright (c) 2003-2009 by University of Queensland
5  # Earth Systems Science Computational Center (ESSCC)  # Earth Systems Science Computational Center (ESSCC)
6  # http://www.uq.edu.au/esscc  # http://www.uq.edu.au/esscc
7  #  #
# Line 11  Line 11 
11  #  #
12  ########################################################  ########################################################
13    
14  __copyright__="""Copyright (c) 2003-2008 by University of Queensland  __copyright__="""Copyright (c) 2003-2009 by University of Queensland
15  Earth Systems Science Computational Center (ESSCC)  Earth Systems Science Computational Center (ESSCC)
16  http://www.uq.edu.au/esscc  http://www.uq.edu.au/esscc
17  Primary Business: Queensland, Australia"""  Primary Business: Queensland, Australia"""
# Line 94  class SolverOptions(object): Line 94  class SolverOptions(object):
94      @cvar YAIR_SHAPIRA_COARSENING: AMG coarsening method by Yair-Shapira      @cvar YAIR_SHAPIRA_COARSENING: AMG coarsening method by Yair-Shapira
95      @cvar RUGE_STUEBEN_COARSENING: AMG coarsening method by Ruge and Stueben      @cvar RUGE_STUEBEN_COARSENING: AMG coarsening method by Ruge and Stueben
96      @cvar AGGREGATION_COARSENING: AMG coarsening using (symmetric) aggregation      @cvar AGGREGATION_COARSENING: AMG coarsening using (symmetric) aggregation
97        @cvar MIN_COARSE_MATRIX_SIZE: minimum size of the coarsest level matrix to use direct solver.
98      @cvar NO_PRECONDITIONER: no preconditioner is applied.      @cvar NO_PRECONDITIONER: no preconditioner is applied.
99      """      """
100      DEFAULT= 0      DEFAULT= 0
# Line 132  class SolverOptions(object): Line 133  class SolverOptions(object):
133      RUGE_STUEBEN_COARSENING=34      RUGE_STUEBEN_COARSENING=34
134      AGGREGATION_COARSENING=35      AGGREGATION_COARSENING=35
135      NO_PRECONDITIONER=36      NO_PRECONDITIONER=36
136        MIN_COARSE_MATRIX_SIZE=37
137        
138      def __init__(self):      def __init__(self):
139          self.setLevelMax()          self.setLevelMax()
140          self.setCoarseningThreshold()          self.setCoarseningThreshold()
# Line 156  class SolverOptions(object): Line 159  class SolverOptions(object):
159          self.setSolverMethod()          self.setSolverMethod()
160          self.setPreconditioner()          self.setPreconditioner()
161          self.setCoarsening()          self.setCoarsening()
162            self.setMinCoarseMatrixSize()
163          self.setRelaxationFactor()                  self.setRelaxationFactor()        
164          self.resetDiagnostics(all=True)          self.resetDiagnostics(all=True)
165    
# Line 172  class SolverOptions(object): Line 176  class SolverOptions(object):
176          out+="\nAbsolute tolerance = %e"%self.getAbsoluteTolerance()          out+="\nAbsolute tolerance = %e"%self.getAbsoluteTolerance()
177          out+="\nSymmetric problem = %s"%self.isSymmetric()          out+="\nSymmetric problem = %s"%self.isSymmetric()
178          out+="\nMaximum number of iteration steps = %s"%self.getIterMax()          out+="\nMaximum number of iteration steps = %s"%self.getIterMax()
179            out+="\nInner tolerance = %e"%self.getInnerTolerance()
180            out+="\nAdapt innner tolerance = %s"%self.adaptInnerTolerance()
181        
182          if self.getPackage() == self.PASO:          if self.getPackage() == self.PASO:
183              out+="\nSolver method = %s"%self.getName(self.getSolverMethod())              out+="\nSolver method = %s"%self.getName(self.getSolverMethod())
184              if self.getSolverMethod() == self.GMRES:              if self.getSolverMethod() == self.GMRES:
# Line 186  class SolverOptions(object): Line 193  class SolverOptions(object):
193              if self.getPreconditioner() == self.AMG:              if self.getPreconditioner() == self.AMG:
194                  out+="\nMaximum number of levels = %s"%self.LevelMax()                  out+="\nMaximum number of levels = %s"%self.LevelMax()
195                  out+="\nCoarsening method = %s"%self.getName(self.getCoarsening())                  out+="\nCoarsening method = %s"%self.getName(self.getCoarsening())
196                  out+="\nCoarsening threshold = %e"%self.getCoarseningThreshold()                  out+="\nCoarsening threshold = %e"%self.getMinCoarseMatrixSize()
197                    out+="\nMinimum size of the coarsest level matrix = %e"%self.getCoarseningThreshold()
198                  out+="\nNumber of pre / post sweeps = %s / %s, %s"%(self.getNumPreSweeps(), self.getNumPostSweeps(), self.getNumSweeps())                  out+="\nNumber of pre / post sweeps = %s / %s, %s"%(self.getNumPreSweeps(), self.getNumPostSweeps(), self.getNumSweeps())
199              if self.getPreconditioner() == self.GAUSS_SEIDEL:              if self.getPreconditioner() == self.GAUSS_SEIDEL:
200                  out+="\nNumber of sweeps = %s"%self.getNumSweeps()                  out+="\nNumber of sweeps = %s"%self.getNumSweeps()
# Line 239  class SolverOptions(object): Line 247  class SolverOptions(object):
247          if key == self.RUGE_STUEBEN_COARSENING: return "RUGE_STUEBEN_COARSENING"          if key == self.RUGE_STUEBEN_COARSENING: return "RUGE_STUEBEN_COARSENING"
248          if key == self.AGGREGATION_COARSENING: return "AGGREGATION_COARSENING"          if key == self.AGGREGATION_COARSENING: return "AGGREGATION_COARSENING"
249          if key == self.NO_PRECONDITIONER: return "NO_PRECONDITIONER"          if key == self.NO_PRECONDITIONER: return "NO_PRECONDITIONER"
250            if key == self.MIN_COARSE_MATRIX_SIZE: return "MIN_COARSE_MATRIX_SIZE"
251                    
252      def resetDiagnostics(self,all=False):      def resetDiagnostics(self,all=False):
253          """          """
# Line 331  class SolverOptions(object): Line 340  class SolverOptions(object):
340          Sets the key of the coarsening method to be applied in AMG.          Sets the key of the coarsening method to be applied in AMG.
341    
342          @param method: selects the coarsening method .          @param method: selects the coarsening method .
343          @type method: in L{SolverOptions.DEFAULT}, L{SolverOptions.YAIR_SHAPIRA_COARSENING},          @type method: in {SolverOptions.DEFAULT}, L{SolverOptions.YAIR_SHAPIRA_COARSENING},
344          L{SolverOptions.RUGE_STUEBEN_COARSENING}, L{SolverOptions.AGGREGATION_COARSENING}          L{SolverOptions.RUGE_STUEBEN_COARSENING}, L{SolverOptions.AGGREGATION_COARSENING}
345          """          """
346        if method==None: method=0
347          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]:
348               raise ValueError,"unknown coarsening method %s"%method               raise ValueError,"unknown coarsening method %s"%method
349          self.__coarsening=method          self.__coarsening=method
350        
351      def getCoarsening(self):      def getCoarsening(self):
352          """          """
353          Returns the key of the coarsening algorithm to be applied AMG.          Returns the key of the coarsening algorithm to be applied AMG.
# Line 345  class SolverOptions(object): Line 356  class SolverOptions(object):
356          L{SolverOptions.RUGE_STUEBEN_COARSENING}, L{SolverOptions.AGGREGATION_COARSENING}          L{SolverOptions.RUGE_STUEBEN_COARSENING}, L{SolverOptions.AGGREGATION_COARSENING}
357          """          """
358          return self.__coarsening          return self.__coarsening
359          
360        def setMinCoarseMatrixSize(self,size=500):
361            """
362            Sets the minumum size of the coarsest level matrix in AMG.
363    
364            @param size: minumum size of the coarsest level matrix .
365            @type size: positive C{int} or C{None}
366            """
367            size=int(size)
368            if size<0:
369               raise ValueError,"minumum size of the coarsest level matrix must be non-negative."
370        if size==None: size=500
371            self.__MinCoarseMatrixSize=size
372            
373        def getMinCoarseMatrixSize(self):
374            """
375            Returns the minumum size of the coarsest level matrix in AMG.
376    
377            @rtype: C{int}
378            """
379            return self.__MinCoarseMatrixSize
380          
381      def setPreconditioner(self, preconditioner=10):      def setPreconditioner(self, preconditioner=10):
382          """          """
383          Sets the preconditioner to be used.          Sets the preconditioner to be used.
# Line 356  class SolverOptions(object): Line 389  class SolverOptions(object):
389          @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
390          an unknown preconditioner.          an unknown preconditioner.
391          """          """
392        if preconditioner==None: preconditioner=10
393          if not preconditioner in [ SolverOptions.SSOR, SolverOptions.ILU0, SolverOptions.ILUT, SolverOptions.JACOBI,          if not preconditioner in [ SolverOptions.SSOR, SolverOptions.ILU0, SolverOptions.ILUT, SolverOptions.JACOBI,
394                                      SolverOptions.AMG, SolverOptions.REC_ILU, SolverOptions.GAUSS_SEIDEL, SolverOptions.RILU,                                      SolverOptions.AMG, SolverOptions.REC_ILU, SolverOptions.GAUSS_SEIDEL, SolverOptions.RILU,
395                                      SolverOptions.NO_PRECONDITIONER] :                                      SolverOptions.NO_PRECONDITIONER] :
# Line 385  class SolverOptions(object): Line 419  class SolverOptions(object):
419          @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
420          an unknown solver method.          an unknown solver method.
421          """          """
422        if method==None: method=0
423          if not method in [ SolverOptions.DEFAULT, SolverOptions.DIRECT, SolverOptions.CHOLEVSKY, SolverOptions.PCG,          if not method in [ SolverOptions.DEFAULT, SolverOptions.DIRECT, SolverOptions.CHOLEVSKY, SolverOptions.PCG,
424                             SolverOptions.CR, SolverOptions.CGS, SolverOptions.BICGSTAB, SolverOptions.SSOR,                             SolverOptions.CR, SolverOptions.CGS, SolverOptions.BICGSTAB, SolverOptions.SSOR,
425                             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 446  class SolverOptions(object):
446          @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}
447          @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.
448          """          """
449        if package==None: package=0
450          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]:
451               raise ValueError,"unknown solver package %s"%package               raise ValueError,"unknown solver package %s"%package
452          self.__package=package          self.__package=package
# Line 456  class SolverOptions(object): Line 492  class SolverOptions(object):
492              if restart<1:              if restart<1:
493                  raise ValueError,"restart must be positive."                  raise ValueError,"restart must be positive."
494              self.__restart=restart              self.__restart=restart
495            
496      def getRestart(self):      def getRestart(self):
497          """          """
498          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 504  class SolverOptions(object):
504              return None              return None
505          else:          else:
506              return self.__restart              return self.__restart
507        def _getRestartForC(self):
508            r=self.getRestart()
509            if r == None:
510                return -1
511                else:
512                return r
513      def setTruncation(self,truncation=20):      def setTruncation(self,truncation=20):
514          """          """
515          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 504  class SolverOptions(object): Line 547  class SolverOptions(object):
547          @rtype: C{int}          @rtype: C{int}
548          """          """
549          return self.__inner_iter_max          return self.__inner_iter_max
550      def setIterMax(self,iter_max=10000):      def setIterMax(self,iter_max=100000):
551          """          """
552          Sets the maximum number of iteration steps          Sets the maximum number of iteration steps
553    
# Line 775  class SolverOptions(object): Line 818  class SolverOptions(object):
818          Switches the verbosity of the solver off.          Switches the verbosity of the solver off.
819          """          """
820          self.__verbose=False          self.__verbose=False
821      def setVerbosity(self,verbose=True):      def setVerbosity(self,verbose=False):
822          """          """
823          Sets the verbosity flag for the solver to C{flag}.          Sets the verbosity flag for the solver to C{flag}.
824    
# Line 786  class SolverOptions(object): Line 829  class SolverOptions(object):
829              self.setVerbosityOn()              self.setVerbosityOn()
830          else:          else:
831              self.setVerbosityOff()              self.setVerbosityOff()
832            
         self.__adapt_inner_tolerance=True  
833      def adaptInnerTolerance(self):      def adaptInnerTolerance(self):
834          """          """
835          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 1250  class LinearProblem(object):
1250     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
1251     problem will be solved to get the unknown M{u}.     problem will be solved to get the unknown M{u}.
1252    
1253     @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"  
   
   
1254     def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):     def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):
1255       """       """
1256       Initializes a linear problem.       Initializes a linear problem.
# Line 1305  class LinearProblem(object): Line 1274  class LinearProblem(object):
1274       self.__altered_coefficients=False       self.__altered_coefficients=False
1275       self.__reduce_equation_order=False       self.__reduce_equation_order=False
1276       self.__reduce_solution_order=False       self.__reduce_solution_order=False
1277       self.__tolerance=1.e-8       self.__sym=False
1278       self.__solver_method=self.DEFAULT       self.setSolverOptions()
      self.__solver_package=self.DEFAULT  
      self.__preconditioner=self.DEFAULT  
1279       self.__is_RHS_valid=False       self.__is_RHS_valid=False
1280       self.__is_operator_valid=False       self.__is_operator_valid=False
      self.__sym=False  
1281       self.__COEFFICIENTS={}       self.__COEFFICIENTS={}
1282         self.__solution_rtol=1.e99
1283         self.__solution_atol=1.e99
1284         self.setSymmetryOff()
1285       # initialize things:       # initialize things:
1286       self.resetAllCoefficients()       self.resetAllCoefficients()
1287       self.__system_type=None       self.initializeSystem()
      self.updateSystemType()  
1288     # ==========================================================================     # ==========================================================================
1289     #    general stuff:     #    general stuff:
1290     # ==========================================================================     # ==========================================================================
# Line 1392  class LinearProblem(object): Line 1360  class LinearProblem(object):
1360       @rtype: L{Domain<escript.Domain>}       @rtype: L{Domain<escript.Domain>}
1361       """       """
1362       return self.__domain       return self.__domain
1363       def getDomainStatus(self):
1364         """
1365         Return the status indicator of the domain
1366         """
1367         return self.getDomain().getStatus()
1368    
1369       def getSystemStatus(self):
1370         """
1371         Return the domain status used to build the current system
1372         """
1373         return self.__system_status
1374       def setSystemStatus(self,status=None):
1375         """
1376         Sets the system status to C{status} if C{status} is not present the
1377         current status of the domain is used.
1378         """
1379         if status == None:
1380             self.__system_status=self.getDomainStatus()
1381         else:
1382             self.__system_status=status
1383    
1384     def getDim(self):     def getDim(self):
1385       """       """
# Line 1481  class LinearProblem(object): Line 1469  class LinearProblem(object):
1469     # ==========================================================================     # ==========================================================================
1470     #   solver settings:     #   solver settings:
1471     # ==========================================================================     # ==========================================================================
1472     def setSolverMethod(self,solver=None,preconditioner=None):     def setSolverOptions(self,options=None):
1473         """         """
1474         Sets a new solver method and/or preconditioner.         Sets the solver options.
   
        @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}  
        """  
        return self.__solver_package  
1475    
1476           @param options: the new solver options. If equal C{None}, the solver options are set to the default.
1477           @type options: L{SolverOptions} or C{None}
1478           @note: The symmetry flag of options is overwritten by the symmetry flag of the L{LinearProblem}.
1479           """
1480           if options==None:
1481              self.__solver_options=SolverOptions()
1482           elif isinstance(options, SolverOptions):
1483              self.__solver_options=options
1484           else:
1485              raise ValueError,"options must be a SolverOptions object."
1486           self.__solver_options.setSymmetry(self.__sym)
1487        
1488       def getSolverOptions(self):
1489           """
1490           Returns the solver options
1491      
1492           @rtype: L{SolverOptions}
1493           """
1494           self.__solver_options.setSymmetry(self.__sym)
1495           return self.__solver_options
1496          
1497     def isUsingLumping(self):     def isUsingLumping(self):
1498        """        """
1499        Checks if matrix lumping is the current solver method.        Checks if matrix lumping is the current solver method.
# Line 1588  class LinearProblem(object): Line 1501  class LinearProblem(object):
1501        @return: True if the current solver method is lumping        @return: True if the current solver method is lumping
1502        @rtype: C{bool}        @rtype: C{bool}
1503        """        """
1504        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  
   
1505     # ==========================================================================     # ==========================================================================
1506     #    symmetry  flag:     #    symmetry  flag:
1507     # ==========================================================================     # ==========================================================================
# Line 1624  class LinearProblem(object): Line 1511  class LinearProblem(object):
1511    
1512        @return: True if a symmetric PDE is indicated, False otherwise        @return: True if a symmetric PDE is indicated, False otherwise
1513        @rtype: C{bool}        @rtype: C{bool}
1514          @note: the method is equivalent to use getSolverOptions().isSymmetric()
1515        """        """
1516        return self.__sym        self.getSolverOptions().isSymmetric()
1517    
1518     def setSymmetryOn(self):     def setSymmetryOn(self):
1519        """        """
1520        Sets the symmetry flag.        Sets the symmetry flag.
1521          @note: The method overwrites the symmetry flag set by the solver options
1522        """        """
1523        if not self.isSymmetric():        self.__sym=True
1524           self.trace("PDE is set to be symmetric")        self.getSolverOptions().setSymmetryOn()
          self.__sym=True  
          self.updateSystemType()  
1525    
1526     def setSymmetryOff(self):     def setSymmetryOff(self):
1527        """        """
1528        Clears the symmetry flag.        Clears the symmetry flag.
1529          @note: The method overwrites the symmetry flag set by the solver options
1530        """        """
1531        if self.isSymmetric():        self.__sym=False
1532           self.trace("PDE is set to be nonsymmetric")        self.getSolverOptions().setSymmetryOff()
          self.__sym=False  
          self.updateSystemType()  
1533    
1534     def setSymmetryTo(self,flag=False):     def setSymmetry(self,flag=False):
1535        """        """
1536        Sets the symmetry flag to C{flag}.        Sets the symmetry flag to C{flag}.
1537    
1538        @param flag: If True, the symmetry flag is set otherwise reset.        @param flag: If True, the symmetry flag is set otherwise reset.
1539        @type flag: C{bool}        @type flag: C{bool}
1540          @note: The method overwrites the symmetry flag set by the solver options
1541        """        """
1542        if flag:        self.getSolverOptions().setSymmetry(flag)
          self.setSymmetryOn()  
       else:  
          self.setSymmetryOff()  
   
1543     # ==========================================================================     # ==========================================================================
1544     # function space handling for the equation as well as the solution     # function space handling for the equation as well as the solution
1545     # ==========================================================================     # ==========================================================================
# Line 1783  class LinearProblem(object): Line 1666  class LinearProblem(object):
1666          self.setReducedOrderForEquationOn()          self.setReducedOrderForEquationOn()
1667       else:       else:
1668          self.setReducedOrderForEquationOff()          self.setReducedOrderForEquationOff()
1669       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):  
1670        """        """
1671        Returns the current system type.        Returns the current system type.
1672        """        """
1673        return self.__system_type        return self.__operator_type
1674    
1675     def checkSymmetricTensor(self,name,verbose=True):     def checkSymmetricTensor(self,name,verbose=True):
1676        """        """
# Line 1813  class LinearProblem(object): Line 1684  class LinearProblem(object):
1684        @return: True if coefficient C{name} is symmetric        @return: True if coefficient C{name} is symmetric
1685        @rtype: C{bool}        @rtype: C{bool}
1686        """        """
1687          SMALL_TOLERANCE=util.EPSILON*10.
1688        A=self.getCoefficient(name)        A=self.getCoefficient(name)
1689        verbose=verbose or self.__debug        verbose=verbose or self.__debug
1690        out=True        out=True
1691        if not A.isEmpty():        if not A.isEmpty():
1692           tol=util.Lsup(A)*self.SMALL_TOLERANCE           tol=util.Lsup(A)*SMALL_TOLERANCE
1693           s=A.getShape()           s=A.getShape()
1694           if A.getRank() == 4:           if A.getRank() == 4:
1695              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 1734  class LinearProblem(object):
1734                 symmetric.                 symmetric.
1735        @rtype: C{bool}        @rtype: C{bool}
1736        """        """
1737          SMALL_TOLERANCE=util.EPSILON*10.
1738        B=self.getCoefficient(name0)        B=self.getCoefficient(name0)
1739        C=self.getCoefficient(name1)        C=self.getCoefficient(name1)
1740        verbose=verbose or self.__debug        verbose=verbose or self.__debug
# Line 1875  class LinearProblem(object): Line 1748  class LinearProblem(object):
1748        elif not B.isEmpty() and not C.isEmpty():        elif not B.isEmpty() and not C.isEmpty():
1749           sB=B.getShape()           sB=B.getShape()
1750           sC=C.getShape()           sC=C.getShape()
1751           tol=(util.Lsup(B)+util.Lsup(C))*self.SMALL_TOLERANCE/2.           tol=(util.Lsup(B)+util.Lsup(C))*SMALL_TOLERANCE/2.
1752           if len(sB) != len(sC):           if len(sB) != len(sC):
1753               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))
1754               out=False               out=False
# Line 2020  class LinearProblem(object): Line 1893  class LinearProblem(object):
1893         """         """
1894         Returns True if the solution is still valid.         Returns True if the solution is still valid.
1895         """         """
1896           if not self.getDomainStatus()==self.getSystemStatus(): self.invalidateSolution()
1897           if self.__solution_rtol>self.getSolverOptions().getTolerance() or \
1898              self.__solution_atol>self.getSolverOptions().getAbsoluteTolerance():
1899             self.invalidateSolution()  
1900         return self.__is_solution_valid         return self.__is_solution_valid
1901    
1902     def validOperator(self):     def validOperator(self):
# Line 2040  class LinearProblem(object): Line 1917  class LinearProblem(object):
1917         """         """
1918         Returns True if the operator is still valid.         Returns True if the operator is still valid.
1919         """         """
1920           if not self.getDomainStatus()==self.getSystemStatus(): self.invalidateOperator()
1921           if not self.getRequiredOperatorType()==self.getOperatorType(): self.invalidateOperator()
1922         return self.__is_operator_valid         return self.__is_operator_valid
1923    
1924     def validRightHandSide(self):     def validRightHandSide(self):
# Line 2052  class LinearProblem(object): Line 1931  class LinearProblem(object):
1931         """         """
1932         Indicates the right hand side has to be rebuilt next time it is used.         Indicates the right hand side has to be rebuilt next time it is used.
1933         """         """
1934         if self.isRightHandSideValid(): self.trace("Right hand side has to be rebuilt.")         self.trace("Right hand side has to be rebuilt.")
1935         self.invalidateSolution()         self.invalidateSolution()
1936         self.__is_RHS_valid=False         self.__is_RHS_valid=False
1937    
# Line 2060  class LinearProblem(object): Line 1939  class LinearProblem(object):
1939         """         """
1940         Returns True if the operator is still valid.         Returns True if the operator is still valid.
1941         """         """
1942           if not self.getDomainStatus()==self.getSystemStatus(): self.invalidateRightHandSide()
1943         return self.__is_RHS_valid         return self.__is_RHS_valid
1944    
1945     def invalidateSystem(self):     def invalidateSystem(self):
1946         """         """
1947         Announces that everything has to be rebuilt.         Announces that everything has to be rebuilt.
1948         """         """
        if self.isRightHandSideValid(): self.trace("System has to be rebuilt.")  
1949         self.invalidateSolution()         self.invalidateSolution()
1950         self.invalidateOperator()         self.invalidateOperator()
1951         self.invalidateRightHandSide()         self.invalidateRightHandSide()
# Line 2082  class LinearProblem(object): Line 1961  class LinearProblem(object):
1961         Resets the system clearing the operator, right hand side and solution.         Resets the system clearing the operator, right hand side and solution.
1962         """         """
1963         self.trace("New System has been created.")         self.trace("New System has been created.")
1964           self.__operator_type=None
1965           self.setSystemStatus()
1966         self.__operator=escript.Operator()         self.__operator=escript.Operator()
1967         self.__righthandside=escript.Data()         self.__righthandside=escript.Data()
1968         self.__solution=escript.Data()         self.__solution=escript.Data()
# Line 2136  class LinearProblem(object): Line 2017  class LinearProblem(object):
2017    
2018     def setSolution(self,u):     def setSolution(self,u):
2019         """         """
2020         Sets the solution assuming that makes the system valid.         Sets the solution assuming that makes the system valid with the tolrance
2021           defined by the solver options
2022         """         """
2023           self.__solution_rtol=self.getSolverOptions().getTolerance()
2024           self.__solution_atol=self.getSolverOptions().getAbsoluteTolerance()
2025         self.__solution=u         self.__solution=u
2026         self.validSolution()         self.validSolution()
2027    
# Line 2170  class LinearProblem(object): Line 2054  class LinearProblem(object):
2054         Makes sure that the operator is instantiated and returns it initialized         Makes sure that the operator is instantiated and returns it initialized
2055         with zeros.         with zeros.
2056         """         """
2057         if self.__operator.isEmpty():         if self.getOperatorType() == None:
2058             if self.isUsingLumping():             if self.isUsingLumping():
2059                 self.__operator=self.createSolution()                 self.__operator=self.createSolution()
2060             else:             else:
2061                 self.__operator=self.createOperator()                 self.__operator=self.createOperator()
2062           self.__operator_type=self.getRequiredOperatorType()
2063         else:         else:
2064             if self.isUsingLumping():             if self.isUsingLumping():
2065                 self.__operator.setToZero()                 self.__operator.setToZero()
# Line 2186  class LinearProblem(object): Line 2071  class LinearProblem(object):
2071         """         """
2072         Returns the operator in its current state.         Returns the operator in its current state.
2073         """         """
        if self.__operator.isEmpty(): self.resetOperator()  
2074         return self.__operator         return self.__operator
2075    
2076     def setValue(self,**coefficients):     def setValue(self,**coefficients):
# Line 2248  class LinearProblem(object): Line 2132  class LinearProblem(object):
2132     # methods that are typically overwritten when implementing a particular     # methods that are typically overwritten when implementing a particular
2133     # linear problem     # linear problem
2134     # ==========================================================================     # ==========================================================================
2135     def getRequiredSystemType(self):     def getRequiredOperatorType(self):
2136        """        """
2137        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.
2138    
2139        @note: Typically this method is overwritten when implementing a        @note: Typically this method is overwritten when implementing a
2140               particular linear problem.               particular linear problem.
2141        """        """
2142        return 0        return None
2143    
2144     def createOperator(self):     def createOperator(self):
2145         """         """
# Line 2304  class LinearProblem(object): Line 2188  class LinearProblem(object):
2188                linear problem.                linear problem.
2189         """         """
2190         return (self.getCurrentOperator(), self.getCurrentRightHandSide())         return (self.getCurrentOperator(), self.getCurrentRightHandSide())
 #=====================  
2191    
2192  class LinearPDE(LinearProblem):  class LinearPDE(LinearProblem):
2193     """     """
# Line 2489  class LinearPDE(LinearProblem): Line 2372  class LinearPDE(LinearProblem):
2372       """       """
2373       return "<LinearPDE %d>"%id(self)       return "<LinearPDE %d>"%id(self)
2374    
2375     def getRequiredSystemType(self):     def getRequiredOperatorType(self):
2376        """        """
2377        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.
2378        """        """
2379        return self.getDomain().getSystemMatrixTypeId(self.getSolverMethod()[0],self.getSolverMethod()[1],self.getSolverPackage(),self.isSymmetric())        solver_options=self.getSolverOptions()
2380          return self.getDomain().getSystemMatrixTypeId(solver_options.getSolverMethod(), solver_options.getPreconditioner(),solver_options.getPackage(), solver_options.isSymmetric())
2381    
2382     def checkSymmetry(self,verbose=True):     def checkSymmetry(self,verbose=True):
2383        """        """
# Line 2524  class LinearPDE(LinearProblem): Line 2408  class LinearPDE(LinearProblem):
2408         """         """
2409         Returns an instance of a new operator.         Returns an instance of a new operator.
2410         """         """
2411         self.trace("New operator is allocated.")         optype=self.getRequiredOperatorType()
2412           self.trace("New operator of type %s is allocated."%optype)
2413         return self.getDomain().newOperator( \         return self.getDomain().newOperator( \
2414                             self.getNumEquations(), \                             self.getNumEquations(), \
2415                             self.getFunctionSpaceForEquation(), \                             self.getFunctionSpaceForEquation(), \
2416                             self.getNumSolutions(), \                             self.getNumSolutions(), \
2417                             self.getFunctionSpaceForSolution(), \                             self.getFunctionSpaceForSolution(), \
2418                             self.getSystemType())                             optype)
2419    
2420     def getSolution(self,**options):     def getSolution(self):
2421         """         """
2422         Returns the solution of the PDE.         Returns the solution of the PDE.
2423    
2424         @return: the solution         @return: the solution
2425         @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}  
2426         """         """
2427           option_class=self.getSolverOptions()
2428         if not self.isSolutionValid():         if not self.isSolutionValid():
2429            mat,f=self.getSystem()            mat,f=self.getSystem()
2430            if self.isUsingLumping():            if self.isUsingLumping():
2431               self.setSolution(f*1/mat)               self.setSolution(f*1/mat)
2432            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()  
2433               self.trace("PDE is resolved.")               self.trace("PDE is resolved.")
2434               self.trace("solver options: %s"%str(options))               self.trace("solver options: %s"%str(option_class))
2435               self.setSolution(mat.solve(f,options))               self.setSolution(mat.solve(f,option_class))
2436         return self.getCurrentSolution()         return self.getCurrentSolution()
2437    
2438     def getSystem(self):     def getSystem(self):
# Line 2738  class LinearPDE(LinearProblem): Line 2608  class LinearPDE(LinearProblem):
2608                   self.insertConstraint(rhs_only=False)                   self.insertConstraint(rhs_only=False)
2609                   self.trace("New operator has been built.")                   self.trace("New operator has been built.")
2610                   self.validOperator()                   self.validOperator()
2611           self.setSystemStatus()
2612           self.trace("System status is %s."%self.getSystemStatus())
2613         return (self.getCurrentOperator(), self.getCurrentRightHandSide())         return (self.getCurrentOperator(), self.getCurrentRightHandSide())
2614    
2615     def insertConstraint(self, rhs_only=False):     def insertConstraint(self, rhs_only=False):
# Line 3535  class TransportPDE(LinearProblem): Line 3407  class TransportPDE(LinearProblem):
3407           theta=1.           theta=1.
3408         else:         else:
3409           theta=0.5           theta=0.5
3410           optype=self.getRequiredOperatorType()
3411         self.trace("New Transport problem is allocated.")         self.trace("New Transport problem pf type %s is allocated."%optype)
3412         return self.getDomain().newTransportProblem( \         return self.getDomain().newTransportProblem( \
3413                                 theta,                                 theta,
3414                                 self.getNumEquations(), \                                 self.getNumEquations(), \
3415                                 self.getFunctionSpaceForSolution(), \                                 self.getFunctionSpaceForSolution(), \
3416                                 self.getSystemType())                                 optype)
3417    
3418     def setInitialSolution(self,u):     def setInitialSolution(self,u):
3419         """         """
# Line 3561  class TransportPDE(LinearProblem): Line 3433  class TransportPDE(LinearProblem):
3433                raise ValueError,"Illegal shape %s of initial solution."%(u2.getShape(),)                raise ValueError,"Illegal shape %s of initial solution."%(u2.getShape(),)
3434         self.getOperator().setInitialValue(u2)         self.getOperator().setInitialValue(u2)
3435    
3436     def getRequiredSystemType(self):     def getRequiredOperatorType(self):
3437        """        """
3438        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.
3439    
3440        @return: a code to indicate the type of transport problem scheme used        @return: a code to indicate the type of transport problem scheme used
3441        @rtype: C{float}        @rtype: C{float}
3442        """        """
3443        return self.getDomain().getTransportTypeId(self.getSolverMethod()[0],self.getSolverMethod()[1],self.getSolverPackage(),self.isSymmetric())        solver_options=self.getSolverOptions()
3444          return self.getDomain().getTransportTypeId(solver_options.getSolverMethod(), solver_options.getPreconditioner(),solver_options.getPackage(), solver_options.isSymmetric())
3445    
3446     def getUnlimitedTimeStepSize(self):     def getUnlimitedTimeStepSize(self):
3447        """        """
# Line 3613  class TransportPDE(LinearProblem): Line 3486  class TransportPDE(LinearProblem):
3486         """         """
3487         return self.__constraint_factor         return self.__constraint_factor
3488     #====================================================================     #====================================================================
3489     def getSolution(self,dt,**options):     def getSolution(self,dt):
3490         """         """
3491         Returns the solution of the problem.         Returns the solution of the problem.
3492    
3493         @return: the solution         @return: the solution
3494         @rtype: L{Data<escript.Data>}         @rtype: L{Data<escript.Data>}
   
3495         """         """
3496           option_class=self.getSolverOptions()
3497         if dt<=0:         if dt<=0:
3498             raise ValueError,"step size needs to be positive."             raise ValueError,"step size needs to be positive."
3499         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))  
3500         self.validSolution()         self.validSolution()
3501         return self.getCurrentSolution()         return self.getCurrentSolution()
3502    
# Line 3674  class TransportPDE(LinearProblem): Line 3546  class TransportPDE(LinearProblem):
3546            self.trace("New system has been built.")            self.trace("New system has been built.")
3547            self.validOperator()            self.validOperator()
3548            self.validRightHandSide()            self.validRightHandSide()
3549           self.setSystemStatus()
3550           self.trace("System status is %s."%self.getSystemStatus())
3551         return (self.getCurrentOperator(), self.getCurrentRightHandSide())         return (self.getCurrentOperator(), self.getCurrentRightHandSide())
3552    
3553     def setDebug(self, flag):     def setDebug(self, flag):

Legend:
Removed from v.2470  
changed lines
  Added in v.2549

  ViewVC Help
Powered by ViewVC 1.1.26