/[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 3788 by caltinay, Wed Sep 21 04:48:06 2011 UTC revision 3789 by caltinay, Tue Jan 31 04:55:05 2012 UTC
# Line 38  by its advective terms. Line 38  by its advective terms.
38  """  """
39    
40  import math  import math
41  import escript  from . import escript
42  import util  from . import util
43  import numpy  import numpy
44    
45  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
# Line 421  class SolverOptions(object): Line 421  class SolverOptions(object):
421          elif name == "coarse_level_sparsity": return self.__coarse_level_sparsity          elif name == "coarse_level_sparsity": return self.__coarse_level_sparsity
422          elif name == "num_coarse_unknowns": return self.__num_coarse_unknowns          elif name == "num_coarse_unknowns": return self.__num_coarse_unknowns
423          else:          else:
424              raise ValueError,"unknown diagnostic item %s"%name              raise ValueError("unknown diagnostic item %s"%name)
425      def hasConverged(self):      def hasConverged(self):
426          """          """
427          Returns ``True`` if the last solver call has been finalized successfully.          Returns ``True`` if the last solver call has been finalized successfully.
# Line 437  class SolverOptions(object): Line 437  class SolverOptions(object):
437          """          """
438      if method==None: method=0      if method==None: method=0
439          if not method in [self.DEFAULT, self.YAIR_SHAPIRA_COARSENING, self.RUGE_STUEBEN_COARSENING, self.AGGREGATION_COARSENING, self.STANDARD_COARSENING, self.CIJP_FIXED_RANDOM_COARSENING, self.CIJP_COARSENING, self.FALGOUT_COARSENING, self.PMIS_COARSENING, self.HMIS_COARSENING,]:          if not method in [self.DEFAULT, self.YAIR_SHAPIRA_COARSENING, self.RUGE_STUEBEN_COARSENING, self.AGGREGATION_COARSENING, self.STANDARD_COARSENING, self.CIJP_FIXED_RANDOM_COARSENING, self.CIJP_COARSENING, self.FALGOUT_COARSENING, self.PMIS_COARSENING, self.HMIS_COARSENING,]:
440               raise ValueError,"unknown coarsening method %s"%method               raise ValueError("unknown coarsening method %s"%method)
441          self.__coarsening=method          self.__coarsening=method
442            
443      def getCoarsening(self):      def getCoarsening(self):
# Line 458  class SolverOptions(object): Line 458  class SolverOptions(object):
458      if size==None: size=500      if size==None: size=500
459          size=int(size)          size=int(size)
460          if size<0:          if size<0:
461             raise ValueError,"minumum size of the coarsest level matrix must be non-negative."             raise ValueError("minumum size of the coarsest level matrix must be non-negative.")
462          self.__MinCoarseMatrixSize=size          self.__MinCoarseMatrixSize=size
463                    
464      def getMinCoarseMatrixSize(self):      def getMinCoarseMatrixSize(self):
# Line 483  class SolverOptions(object): Line 483  class SolverOptions(object):
483          if not preconditioner in [  SolverOptions.ILU0, SolverOptions.ILUT, SolverOptions.JACOBI,          if not preconditioner in [  SolverOptions.ILU0, SolverOptions.ILUT, SolverOptions.JACOBI,
484                                      SolverOptions.AMG, SolverOptions.AMLI, SolverOptions.REC_ILU, SolverOptions.GAUSS_SEIDEL, SolverOptions.RILU, SolverOptions.BOOMERAMG,                                      SolverOptions.AMG, SolverOptions.AMLI, SolverOptions.REC_ILU, SolverOptions.GAUSS_SEIDEL, SolverOptions.RILU, SolverOptions.BOOMERAMG,
485                                      SolverOptions.NO_PRECONDITIONER] :                                      SolverOptions.NO_PRECONDITIONER] :
486               raise ValueError,"unknown preconditioner %s"%preconditioner               raise ValueError("unknown preconditioner %s"%preconditioner)
487          self.__preconditioner=preconditioner              self.__preconditioner=preconditioner    
488      def getPreconditioner(self):      def getPreconditioner(self):
489          """          """
# Line 504  class SolverOptions(object): Line 504  class SolverOptions(object):
504          """          """
505      if smoother==None: smoother=28      if smoother==None: smoother=28
506          if not smoother in [ SolverOptions.JACOBI, SolverOptions.GAUSS_SEIDEL ] :          if not smoother in [ SolverOptions.JACOBI, SolverOptions.GAUSS_SEIDEL ] :
507               raise ValueError,"unknown smoother %s"%smoother               raise ValueError("unknown smoother %s"%smoother)
508          self.__smoother=smoother              self.__smoother=smoother    
509      def getSmoother(self):      def getSmoother(self):
510          """          """
# Line 533  class SolverOptions(object): Line 533  class SolverOptions(object):
533                             SolverOptions.GMRES, SolverOptions.PRES20, SolverOptions.ROWSUM_LUMPING, SolverOptions.HRZ_LUMPING,                             SolverOptions.GMRES, SolverOptions.PRES20, SolverOptions.ROWSUM_LUMPING, SolverOptions.HRZ_LUMPING,
534                             SolverOptions.ITERATIVE,                             SolverOptions.ITERATIVE,
535                             SolverOptions.NONLINEAR_GMRES, SolverOptions.TFQMR, SolverOptions.MINRES ]:                             SolverOptions.NONLINEAR_GMRES, SolverOptions.TFQMR, SolverOptions.MINRES ]:
536               raise ValueError,"unknown solver method %s"%method               raise ValueError("unknown solver method %s"%method)
537          self.__method=method          self.__method=method
538      def getSolverMethod(self):      def getSolverMethod(self):
539          """          """
# Line 556  class SolverOptions(object): Line 556  class SolverOptions(object):
556          """          """
557      if package==None: package=0      if package==None: package=0
558          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]:
559               raise ValueError,"unknown solver package %s"%package               raise ValueError("unknown solver package %s"%package)
560          self.__package=package          self.__package=package
561      def getPackage(self):      def getPackage(self):
562          """          """
# Line 574  class SolverOptions(object): Line 574  class SolverOptions(object):
574          :type ordering: in 'SolverOptions.NO_REORDERING', 'SolverOptions.MINIMUM_FILL_IN', 'SolverOptions.NESTED_DISSECTION', 'SolverOptions.DEFAULT_REORDERING'          :type ordering: in 'SolverOptions.NO_REORDERING', 'SolverOptions.MINIMUM_FILL_IN', 'SolverOptions.NESTED_DISSECTION', 'SolverOptions.DEFAULT_REORDERING'
575          """          """
576          if not ordering in [self.NO_REORDERING, self.MINIMUM_FILL_IN, self.NESTED_DISSECTION, self.DEFAULT_REORDERING]:          if not ordering in [self.NO_REORDERING, self.MINIMUM_FILL_IN, self.NESTED_DISSECTION, self.DEFAULT_REORDERING]:
577               raise ValueError,"unknown reordering strategy %s"%ordering               raise ValueError("unknown reordering strategy %s"%ordering)
578          self.__reordering=ordering          self.__reordering=ordering
579      def getReordering(self):      def getReordering(self):
580          """          """
# Line 595  class SolverOptions(object): Line 595  class SolverOptions(object):
595          else:          else:
596              restart=int(restart)              restart=int(restart)
597              if restart<1:              if restart<1:
598                  raise ValueError,"restart must be positive."                  raise ValueError("restart must be positive.")
599              self.__restart=restart              self.__restart=restart
600                    
601      def getRestart(self):      def getRestart(self):
# Line 625  class SolverOptions(object): Line 625  class SolverOptions(object):
625          """          """
626          value=float(value)          value=float(value)
627          if value<0 or value>1.:          if value<0 or value>1.:
628         raise ValueError,"Diagonal dominance threshold must be between 0 and 1."         raise ValueError("Diagonal dominance threshold must be between 0 and 1.")
629      self.__diagonal_dominance_threshold=value      self.__diagonal_dominance_threshold=value
630                    
631      def getDiagonalDominanceThreshold(self):      def getDiagonalDominanceThreshold(self):
# Line 646  class SolverOptions(object): Line 646  class SolverOptions(object):
646          """          """
647          truncation=int(truncation)          truncation=int(truncation)
648          if truncation<1:          if truncation<1:
649             raise ValueError,"truncation must be positive."             raise ValueError("truncation must be positive.")
650          self.__truncation=truncation          self.__truncation=truncation
651      def getTruncation(self):      def getTruncation(self):
652          """          """
# Line 664  class SolverOptions(object): Line 664  class SolverOptions(object):
664          """          """
665          iter_max=int(iter_max)          iter_max=int(iter_max)
666          if iter_max<1:          if iter_max<1:
667             raise ValueError,"maximum number of inner iteration must be positive."             raise ValueError("maximum number of inner iteration must be positive.")
668          self.__inner_iter_max=iter_max          self.__inner_iter_max=iter_max
669      def getInnerIterMax(self):      def getInnerIterMax(self):
670          """          """
# Line 682  class SolverOptions(object): Line 682  class SolverOptions(object):
682          """          """
683          iter_max=int(iter_max)          iter_max=int(iter_max)
684          if iter_max<1:          if iter_max<1:
685             raise ValueError,"maximum number of iteration steps must be positive."             raise ValueError("maximum number of iteration steps must be positive.")
686          self.__iter_max=iter_max          self.__iter_max=iter_max
687      def getIterMax(self):      def getIterMax(self):
688          """          """
# Line 700  class SolverOptions(object): Line 700  class SolverOptions(object):
700          """          """
701          level_max=int(level_max)          level_max=int(level_max)
702          if level_max<0:          if level_max<0:
703             raise ValueError,"maximum number of coarsening levels must be non-negative."             raise ValueError("maximum number of coarsening levels must be non-negative.")
704          self.__level_max=level_max          self.__level_max=level_max
705      def getLevelMax(self):      def getLevelMax(self):
706          """          """
# Line 734  class SolverOptions(object): Line 734  class SolverOptions(object):
734          """          """
735          theta=float(theta)          theta=float(theta)
736          if theta<0 or theta>1:          if theta<0 or theta>1:
737             raise ValueError,"threshold must be non-negative and less or equal 1."             raise ValueError("threshold must be non-negative and less or equal 1.")
738          self.__coarsening_threshold=theta          self.__coarsening_threshold=theta
739      def getCoarseningThreshold(self):      def getCoarseningThreshold(self):
740          """          """
# Line 752  class SolverOptions(object): Line 752  class SolverOptions(object):
752          """          """
753          sweeps=int(sweeps)          sweeps=int(sweeps)
754          if sweeps<1:          if sweeps<1:
755             raise ValueError,"number of sweeps must be positive."             raise ValueError("number of sweeps must be positive.")
756          self.__sweeps=sweeps          self.__sweeps=sweeps
757      def getNumSweeps(self):      def getNumSweeps(self):
758          """          """
# Line 770  class SolverOptions(object): Line 770  class SolverOptions(object):
770          """          """
771          sweeps=int(sweeps)          sweeps=int(sweeps)
772          if sweeps<1:          if sweeps<1:
773             raise ValueError,"number of sweeps must be positive."             raise ValueError("number of sweeps must be positive.")
774          self.__pre_sweeps=sweeps          self.__pre_sweeps=sweeps
775      def getNumPreSweeps(self):      def getNumPreSweeps(self):
776          """          """
# Line 788  class SolverOptions(object): Line 788  class SolverOptions(object):
788          """          """
789          sweeps=int(sweeps)          sweeps=int(sweeps)
790          if sweeps<1:          if sweeps<1:
791             raise ValueError,"number of sweeps must be positive."             raise ValueError("number of sweeps must be positive.")
792          self.__post_sweeps=sweeps          self.__post_sweeps=sweeps
793      def getNumPostSweeps(self):      def getNumPostSweeps(self):
794          """          """
# Line 807  class SolverOptions(object): Line 807  class SolverOptions(object):
807          """          """
808          rtol=float(rtol)          rtol=float(rtol)
809          if rtol<0 or rtol>1:          if rtol<0 or rtol>1:
810             raise ValueError,"tolerance must be non-negative and less or equal 1."             raise ValueError("tolerance must be non-negative and less or equal 1.")
811          self.__tolerance=rtol          self.__tolerance=rtol
812      def getTolerance(self):      def getTolerance(self):
813          """          """
# Line 825  class SolverOptions(object): Line 825  class SolverOptions(object):
825          """          """
826          atol=float(atol)          atol=float(atol)
827          if atol<0:          if atol<0:
828             raise ValueError,"tolerance must be non-negative."             raise ValueError("tolerance must be non-negative.")
829          self.__absolute_tolerance=atol          self.__absolute_tolerance=atol
830      def getAbsoluteTolerance(self):      def getAbsoluteTolerance(self):
831          """          """
# Line 844  class SolverOptions(object): Line 844  class SolverOptions(object):
844          """          """
845          rtol=float(rtol)          rtol=float(rtol)
846          if rtol<=0 or rtol>1:          if rtol<=0 or rtol>1:
847              raise ValueError,"tolerance must be positive and less or equal 1."              raise ValueError("tolerance must be positive and less or equal 1.")
848          self.__inner_tolerance=rtol          self.__inner_tolerance=rtol
849      def getInnerTolerance(self):      def getInnerTolerance(self):
850          """          """
# Line 862  class SolverOptions(object): Line 862  class SolverOptions(object):
862          """          """
863          drop_tol=float(drop_tol)          drop_tol=float(drop_tol)
864          if drop_tol<=0 or drop_tol>1:          if drop_tol<=0 or drop_tol>1:
865              raise ValueError,"drop tolerance must be positive and less or equal 1."              raise ValueError("drop tolerance must be positive and less or equal 1.")
866          self.__drop_tolerance=drop_tol          self.__drop_tolerance=drop_tol
867      def getDropTolerance(self):      def getDropTolerance(self):
868          """          """
# Line 881  class SolverOptions(object): Line 881  class SolverOptions(object):
881          """          """
882          storage=float(storage)          storage=float(storage)
883          if storage<1:          if storage<1:
884              raise ValueError,"allowed storage increase must be greater or equal to 1."              raise ValueError("allowed storage increase must be greater or equal to 1.")
885          self.__drop_storage=storage          self.__drop_storage=storage
886      def getDropStorage(self):      def getDropStorage(self):
887            
# Line 901  class SolverOptions(object): Line 901  class SolverOptions(object):
901          """          """
902          factor=float(factor)          factor=float(factor)
903          if factor<0:          if factor<0:
904              raise ValueError,"relaxation factor must be non-negative."              raise ValueError("relaxation factor must be non-negative.")
905          self.__relaxation=factor          self.__relaxation=factor
906      def getRelaxationFactor(self):      def getRelaxationFactor(self):
907            
# Line 1084  class SolverOptions(object): Line 1084  class SolverOptions(object):
1084         """         """
1085         sparsity=float(sparsity)         sparsity=float(sparsity)
1086         if sparsity<0:         if sparsity<0:
1087       raise ValueError,"sparsity must be non-negative."       raise ValueError("sparsity must be non-negative.")
1088         if sparsity>1:         if sparsity>1:
1089            raise ValueError,"sparsity must be less than 1."            raise ValueError("sparsity must be less than 1.")
1090         self.__min_sparsity=sparsity         self.__min_sparsity=sparsity
1091    
1092      def getMinCoarseMatrixSparsity(self):      def getMinCoarseMatrixSparsity(self):
# Line 1109  class SolverOptions(object): Line 1109  class SolverOptions(object):
1109        """        """
1110        refinements=int(refinements)        refinements=int(refinements)
1111        if refinements<0:        if refinements<0:
1112       raise ValueError,"number of refinements must be non-negative."       raise ValueError("number of refinements must be non-negative.")
1113        self.__refinements=refinements        self.__refinements=refinements
1114        
1115      def getNumRefinements(self):      def getNumRefinements(self):
# Line 1129  class SolverOptions(object): Line 1129  class SolverOptions(object):
1129        """        """
1130        refinements=int(refinements)        refinements=int(refinements)
1131        if refinements<0:        if refinements<0:
1132       raise ValueError,"number of refinements must be non-negative."       raise ValueError("number of refinements must be non-negative.")
1133        self.__coarse_refinements=refinements        self.__coarse_refinements=refinements
1134        
1135      def getNumCoarseMatrixRefinements(self):      def getNumCoarseMatrixRefinements(self):
# Line 1183  class SolverOptions(object): Line 1183  class SolverOptions(object):
1183          """          """
1184      if method==None: method=self.DIRECT_INTERPOLATION      if method==None: method=self.DIRECT_INTERPOLATION
1185          if not method in [ SolverOptions.CLASSIC_INTERPOLATION_WITH_FF_COUPLING, SolverOptions.CLASSIC_INTERPOLATION, SolverOptions.DIRECT_INTERPOLATION ]:          if not method in [ SolverOptions.CLASSIC_INTERPOLATION_WITH_FF_COUPLING, SolverOptions.CLASSIC_INTERPOLATION, SolverOptions.DIRECT_INTERPOLATION ]:
1186               raise ValueError,"unknown AMG interpolation method %s"%method               raise ValueError("unknown AMG interpolation method %s"%method)
1187          self.__amg_interpolation_method=method          self.__amg_interpolation_method=method
1188    
1189      def getAMGInterpolation(self):      def getAMGInterpolation(self):
# Line 1388  class PDECoef(object): Line 1388  class PDECoef(object):
1388                  try:                  try:
1389                    newValue=escript.Data(newValue,self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder))                    newValue=escript.Data(newValue,self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder))
1390                  except:                  except:
1391                    raise IllegalCoefficientFunctionSpace,"Unable to interpolate coefficient to function space %s"%self.getFunctionSpace(domain)                    raise IllegalCoefficientFunctionSpace("Unable to interpolate coefficient to function space %s"%self.getFunctionSpace(domain))
1392         else:         else:
1393             newValue=escript.Data(newValue,self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder))             newValue=escript.Data(newValue,self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder))
1394         if not newValue.isEmpty():         if not newValue.isEmpty():
1395             if not self.getShape(domain,numEquations,numSolutions)==newValue.getShape():             if not self.getShape(domain,numEquations,numSolutions)==newValue.getShape():
1396                 raise IllegalCoefficientValue,"Expected shape of coefficient is %s but actual shape is %s."%(self.getShape(domain,numEquations,numSolutions),newValue.getShape())                 raise IllegalCoefficientValue("Expected shape of coefficient is %s but actual shape is %s."%(self.getShape(domain,numEquations,numSolutions),newValue.getShape()))
1397         self.value=newValue         self.value=newValue
1398    
1399      def isAlteringOperator(self):      def isAlteringOperator(self):
# Line 1631  class LinearProblem(object): Line 1631  class LinearProblem(object):
1631       :param text: message to be printed       :param text: message to be printed
1632       :type text: ``string``       :type text: ``string``
1633       """       """
1634       if self.__debug: print "%s: %s"%(str(self),text)       if self.__debug: print(("%s: %s"%(str(self),text)))
1635    
1636     # ==========================================================================     # ==========================================================================
1637     # some service functions:     # some service functions:
# Line 1646  class LinearProblem(object): Line 1646  class LinearProblem(object):
1646    
1647         to introduce the coefficients *A* and *B*.         to introduce the coefficients *A* and *B*.
1648         """         """
1649         for name, type in coeff.items():         for name, type in list(coeff.items()):
1650             if not isinstance(type,PDECoef):             if not isinstance(type,PDECoef):
1651                raise ValueError,"coefficient %s has no type."%name                raise ValueError("coefficient %s has no type."%name)
1652             self.__COEFFICIENTS[name]=type             self.__COEFFICIENTS[name]=type
1653             self.__COEFFICIENTS[name].resetValue()             self.__COEFFICIENTS[name].resetValue()
1654             self.trace("coefficient %s has been introduced."%name)             self.trace("coefficient %s has been introduced."%name)
# Line 1701  class LinearProblem(object): Line 1701  class LinearProblem(object):
1701       """       """
1702       if self.__numEquations==None:       if self.__numEquations==None:
1703           if self.__numSolutions==None:           if self.__numSolutions==None:
1704              raise UndefinedPDEError,"Number of equations is undefined. Please specify argument numEquations."              raise UndefinedPDEError("Number of equations is undefined. Please specify argument numEquations.")
1705           else:           else:
1706              self.__numEquations=self.__numSolutions              self.__numEquations=self.__numSolutions
1707       return self.__numEquations       return self.__numEquations
# Line 1716  class LinearProblem(object): Line 1716  class LinearProblem(object):
1716       """       """
1717       if self.__numSolutions==None:       if self.__numSolutions==None:
1718          if self.__numEquations==None:          if self.__numEquations==None:
1719              raise UndefinedPDEError,"Number of solution is undefined. Please specify argument numSolutions."              raise UndefinedPDEError("Number of solution is undefined. Please specify argument numSolutions.")
1720          else:          else:
1721              self.__numSolutions=self.__numEquations              self.__numSolutions=self.__numEquations
1722       return self.__numSolutions       return self.__numSolutions
# Line 1783  class LinearProblem(object): Line 1783  class LinearProblem(object):
1783         elif isinstance(options, SolverOptions):         elif isinstance(options, SolverOptions):
1784            self.__solver_options=options            self.__solver_options=options
1785         else:         else:
1786            raise ValueError,"options must be a SolverOptions object."            raise ValueError("options must be a SolverOptions object.")
1787         self.__solver_options.setSymmetry(self.__sym)         self.__solver_options.setSymmetry(self.__sym)
1788            
1789     def getSolverOptions(self):     def getSolverOptions(self):
# Line 1889  class LinearProblem(object): Line 1889  class LinearProblem(object):
1889       """       """
1890       if not self.__reduce_solution_order:       if not self.__reduce_solution_order:
1891           if self.__altered_coefficients:           if self.__altered_coefficients:
1892                raise RuntimeError,"order cannot be altered after coefficients have been defined."                raise RuntimeError("order cannot be altered after coefficients have been defined.")
1893           self.trace("Reduced order is used for solution representation.")           self.trace("Reduced order is used for solution representation.")
1894           self.__reduce_solution_order=True           self.__reduce_solution_order=True
1895           self.initializeSystem()           self.initializeSystem()
# Line 1903  class LinearProblem(object): Line 1903  class LinearProblem(object):
1903       """       """
1904       if self.__reduce_solution_order:       if self.__reduce_solution_order:
1905           if self.__altered_coefficients:           if self.__altered_coefficients:
1906                raise RuntimeError,"order cannot be altered after coefficients have been defined."                raise RuntimeError("order cannot be altered after coefficients have been defined.")
1907           self.trace("Full order is used to interpolate solution.")           self.trace("Full order is used to interpolate solution.")
1908           self.__reduce_solution_order=False           self.__reduce_solution_order=False
1909           self.initializeSystem()           self.initializeSystem()
# Line 1933  class LinearProblem(object): Line 1933  class LinearProblem(object):
1933       """       """
1934       if not self.__reduce_equation_order:       if not self.__reduce_equation_order:
1935           if self.__altered_coefficients:           if self.__altered_coefficients:
1936                raise RuntimeError,"order cannot be altered after coefficients have been defined."                raise RuntimeError("order cannot be altered after coefficients have been defined.")
1937           self.trace("Reduced order is used for test functions.")           self.trace("Reduced order is used for test functions.")
1938           self.__reduce_equation_order=True           self.__reduce_equation_order=True
1939           self.initializeSystem()           self.initializeSystem()
# Line 1947  class LinearProblem(object): Line 1947  class LinearProblem(object):
1947       """       """
1948       if self.__reduce_equation_order:       if self.__reduce_equation_order:
1949           if self.__altered_coefficients:           if self.__altered_coefficients:
1950                raise RuntimeError,"order cannot be altered after coefficients have been defined."                raise RuntimeError("order cannot be altered after coefficients have been defined.")
1951           self.trace("Full order is used for test functions.")           self.trace("Full order is used for test functions.")
1952           self.__reduce_equation_order=False           self.__reduce_equation_order=False
1953           self.initializeSystem()           self.initializeSystem()
# Line 1999  class LinearProblem(object): Line 1999  class LinearProblem(object):
1999                       for k in range(s[2]):                       for k in range(s[2]):
2000                          for l in range(s[3]):                          for l in range(s[3]):
2001                              if util.Lsup(A[i,j,k,l]-A[k,l,i,j])>tol:                              if util.Lsup(A[i,j,k,l]-A[k,l,i,j])>tol:
2002                                 if verbose: print "non-symmetric problem as %s[%d,%d,%d,%d]!=%s[%d,%d,%d,%d]"%(name,i,j,k,l,name,k,l,i,j)                                 if verbose: print(("non-symmetric problem as %s[%d,%d,%d,%d]!=%s[%d,%d,%d,%d]"%(name,i,j,k,l,name,k,l,i,j)))
2003                                 out=False                                 out=False
2004              else:              else:
2005                   if verbose: print "non-symmetric problem because of inappropriate shape %s of coefficient %s."%(s,name)                   if verbose: print(("non-symmetric problem because of inappropriate shape %s of coefficient %s."%(s,name)))
2006                   out=False                   out=False
2007           elif A.getRank() == 2:           elif A.getRank() == 2:
2008              if s[0]==s[1]:              if s[0]==s[1]:
2009                 for j in range(s[0]):                 for j in range(s[0]):
2010                    for l in range(s[1]):                    for l in range(s[1]):
2011                       if util.Lsup(A[j,l]-A[l,j])>tol:                       if util.Lsup(A[j,l]-A[l,j])>tol:
2012                          if verbose: print "non-symmetric problem because %s[%d,%d]!=%s[%d,%d]"%(name,j,l,name,l,j)                          if verbose: print(("non-symmetric problem because %s[%d,%d]!=%s[%d,%d]"%(name,j,l,name,l,j)))
2013                          out=False                          out=False
2014              else:              else:
2015                   if verbose: print "non-symmetric problem because of inappropriate shape %s of coefficient %s."%(s,name)                   if verbose: print(("non-symmetric problem because of inappropriate shape %s of coefficient %s."%(s,name)))
2016                   out=False                   out=False
2017           elif A.getRank() == 0:           elif A.getRank() == 0:
2018              pass              pass
2019           else:           else:
2020               raise ValueError,"Cannot check rank %s of %s."%(A.getRank(),name)               raise ValueError("Cannot check rank %s of %s."%(A.getRank(),name))
2021        return out        return out
2022    
2023     def checkReciprocalSymmetry(self,name0,name1,verbose=True):     def checkReciprocalSymmetry(self,name0,name1,verbose=True):
# Line 2041  class LinearProblem(object): Line 2041  class LinearProblem(object):
2041        verbose=verbose or self.__debug        verbose=verbose or self.__debug
2042        out=True        out=True
2043        if B.isEmpty() and not C.isEmpty():        if B.isEmpty() and not C.isEmpty():
2044           if verbose: print "non-symmetric problem because %s is not present but %s is"%(name0,name1)           if verbose: print(("non-symmetric problem because %s is not present but %s is"%(name0,name1)))
2045           out=False           out=False
2046        elif not B.isEmpty() and C.isEmpty():        elif not B.isEmpty() and C.isEmpty():
2047           if verbose: print "non-symmetric problem because %s is not present but %s is"%(name0,name1)           if verbose: print(("non-symmetric problem because %s is not present but %s is"%(name0,name1)))
2048           out=False           out=False
2049        elif not B.isEmpty() and not C.isEmpty():        elif not B.isEmpty() and not C.isEmpty():
2050           sB=B.getShape()           sB=B.getShape()
2051           sC=C.getShape()           sC=C.getShape()
2052           tol=(util.Lsup(B)+util.Lsup(C))*SMALL_TOLERANCE/2.           tol=(util.Lsup(B)+util.Lsup(C))*SMALL_TOLERANCE/2.
2053           if len(sB) != len(sC):           if len(sB) != len(sC):
2054               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))))
2055               out=False               out=False
2056           else:           else:
2057               if len(sB)==0:               if len(sB)==0:
2058                 if util.Lsup(B-C)>tol:                 if util.Lsup(B-C)>tol:
2059                    if verbose: print "non-symmetric problem because %s!=%s"%(name0,name1)                    if verbose: print(("non-symmetric problem because %s!=%s"%(name0,name1)))
2060                    out=False                    out=False
2061               elif len(sB)==1:               elif len(sB)==1:
2062                 if sB[0]==sC[0]:                 if sB[0]==sC[0]:
2063                    for j in range(sB[0]):                    for j in range(sB[0]):
2064                       if util.Lsup(B[j]-C[j])>tol:                       if util.Lsup(B[j]-C[j])>tol:
2065                          if verbose: print "non-symmetric PDE because %s[%d]!=%s[%d]"%(name0,j,name1,j)                          if verbose: print(("non-symmetric PDE because %s[%d]!=%s[%d]"%(name0,j,name1,j)))
2066                          out=False                          out=False
2067                 else:                 else:
2068                   if verbose: print "non-symmetric problem because of inappropriate shapes %s and %s of coefficients %s and %s, respectively."%(sB,sC,name0,name1)                   if verbose: print(("non-symmetric problem because of inappropriate shapes %s and %s of coefficients %s and %s, respectively."%(sB,sC,name0,name1)))
2069               elif len(sB)==3:               elif len(sB)==3:
2070                 if sB[0]==sC[1] and sB[1]==sC[2] and sB[2]==sC[0]:                 if sB[0]==sC[1] and sB[1]==sC[2] and sB[2]==sC[0]:
2071                     for i in range(sB[0]):                     for i in range(sB[0]):
2072                        for j in range(sB[1]):                        for j in range(sB[1]):
2073                           for k in range(sB[2]):                           for k in range(sB[2]):
2074                              if util.Lsup(B[i,j,k]-C[k,i,j])>tol:                              if util.Lsup(B[i,j,k]-C[k,i,j])>tol:
2075                                   if verbose: print "non-symmetric problem because %s[%d,%d,%d]!=%s[%d,%d,%d]"%(name0,i,j,k,name1,k,i,j)                                   if verbose: print(("non-symmetric problem because %s[%d,%d,%d]!=%s[%d,%d,%d]"%(name0,i,j,k,name1,k,i,j)))
2076                                   out=False                                   out=False
2077                 else:                 else:
2078                   if verbose: print "non-symmetric problem because of inappropriate shapes %s and %s of coefficients %s and %s, respectively."%(sB,sC,name0,name1)                   if verbose: print(("non-symmetric problem because of inappropriate shapes %s and %s of coefficients %s and %s, respectively."%(sB,sC,name0,name1)))
2079               else:               else:
2080                   raise ValueError,"Cannot check rank %s of %s and %s."%(len(sB),name0,name1)                   raise ValueError("Cannot check rank %s of %s and %s."%(len(sB),name0,name1))
2081        return out        return out
2082    
2083     def getCoefficient(self,name):     def getCoefficient(self,name):
# Line 2093  class LinearProblem(object): Line 2093  class LinearProblem(object):
2093       if self.hasCoefficient(name):       if self.hasCoefficient(name):
2094           return self.__COEFFICIENTS[name].getValue()           return self.__COEFFICIENTS[name].getValue()
2095       else:       else:
2096          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name          raise IllegalCoefficient("illegal coefficient %s requested for general PDE."%name)
2097    
2098     def hasCoefficient(self,name):     def hasCoefficient(self,name):
2099       """       """
# Line 2105  class LinearProblem(object): Line 2105  class LinearProblem(object):
2105                False otherwise                False otherwise
2106       :rtype: ``bool``       :rtype: ``bool``
2107       """       """
2108       return self.__COEFFICIENTS.has_key(name)       return name in self.__COEFFICIENTS
2109    
2110     def createCoefficient(self, name):     def createCoefficient(self, name):
2111       """       """
# Line 2119  class LinearProblem(object): Line 2119  class LinearProblem(object):
2119       if self.hasCoefficient(name):       if self.hasCoefficient(name):
2120          return escript.Data(0.,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name))          return escript.Data(0.,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name))
2121       else:       else:
2122          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name          raise IllegalCoefficient("illegal coefficient %s requested for general PDE."%name)
2123    
2124     def getFunctionSpaceForCoefficient(self,name):     def getFunctionSpaceForCoefficient(self,name):
2125       """       """
# Line 2135  class LinearProblem(object): Line 2135  class LinearProblem(object):
2135       if self.hasCoefficient(name):       if self.hasCoefficient(name):
2136          return self.__COEFFICIENTS[name].getFunctionSpace(self.getDomain())          return self.__COEFFICIENTS[name].getFunctionSpace(self.getDomain())
2137       else:       else:
2138          raise ValueError,"unknown coefficient %s requested"%name          raise ValueError("unknown coefficient %s requested"%name)
2139    
2140     def getShapeOfCoefficient(self,name):     def getShapeOfCoefficient(self,name):
2141       """       """
# Line 2150  class LinearProblem(object): Line 2150  class LinearProblem(object):
2150       if self.hasCoefficient(name):       if self.hasCoefficient(name):
2151          return self.__COEFFICIENTS[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())          return self.__COEFFICIENTS[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())
2152       else:       else:
2153          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name          raise IllegalCoefficient("illegal coefficient %s requested for general PDE."%name)
2154    
2155     def resetAllCoefficients(self):     def resetAllCoefficients(self):
2156       """       """
2157       Resets all coefficients to their default values.       Resets all coefficients to their default values.
2158       """       """
2159       for i in self.__COEFFICIENTS.iterkeys():       for i in list(self.__COEFFICIENTS.keys()):
2160           self.__COEFFICIENTS[i].resetValue()           self.__COEFFICIENTS[i].resetValue()
2161    
2162     def alteredCoefficient(self,name):     def alteredCoefficient(self,name):
# Line 2175  class LinearProblem(object): Line 2175  class LinearProblem(object):
2175             if self.__COEFFICIENTS[name].isAlteringOperator(): self.invalidateOperator()             if self.__COEFFICIENTS[name].isAlteringOperator(): self.invalidateOperator()
2176             if self.__COEFFICIENTS[name].isAlteringRightHandSide(): self.invalidateRightHandSide()             if self.__COEFFICIENTS[name].isAlteringRightHandSide(): self.invalidateRightHandSide()
2177       else:       else:
2178          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name          raise IllegalCoefficient("illegal coefficient %s requested for general PDE."%name)
2179    
2180     def validSolution(self):     def validSolution(self):
2181         """         """
# Line 2385  class LinearProblem(object): Line 2385  class LinearProblem(object):
2385        :raise IllegalCoefficient: if an unknown coefficient keyword is used        :raise IllegalCoefficient: if an unknown coefficient keyword is used
2386        """        """
2387        # check if the coefficients are  legal:        # check if the coefficients are  legal:
2388        for i in coefficients.iterkeys():        for i in list(coefficients.keys()):
2389           if not self.hasCoefficient(i):           if not self.hasCoefficient(i):
2390              raise IllegalCoefficient,"Attempt to set unknown coefficient %s"%i              raise IllegalCoefficient("Attempt to set unknown coefficient %s"%i)
2391        # if the number of unknowns or equations is still unknown we try to estimate them:        # if the number of unknowns or equations is still unknown we try to estimate them:
2392        if self.__numEquations==None or self.__numSolutions==None:        if self.__numEquations==None or self.__numSolutions==None:
2393           for i,d in coefficients.iteritems():           for i,d in list(coefficients.items()):
2394              if hasattr(d,"shape"):              if hasattr(d,"shape"):
2395                  s=d.shape                  s=d.shape
2396              elif hasattr(d,"getShape"):              elif hasattr(d,"getShape"):
# Line 2401  class LinearProblem(object): Line 2401  class LinearProblem(object):
2401                  # get number of equations and number of unknowns:                  # get number of equations and number of unknowns:
2402                  res=self.__COEFFICIENTS[i].estimateNumEquationsAndNumSolutions(self.getDomain(),s)                  res=self.__COEFFICIENTS[i].estimateNumEquationsAndNumSolutions(self.getDomain(),s)
2403                  if res==None:                  if res==None:
2404                      raise IllegalCoefficientValue,"Illegal shape %s of coefficient %s"%(s,i)                      raise IllegalCoefficientValue("Illegal shape %s of coefficient %s"%(s,i))
2405                  else:                  else:
2406                      if self.__numEquations==None: self.__numEquations=res[0]                      if self.__numEquations==None: self.__numEquations=res[0]
2407                      if self.__numSolutions==None: self.__numSolutions=res[1]                      if self.__numSolutions==None: self.__numSolutions=res[1]
2408        if self.__numEquations==None: raise UndefinedPDEError,"unidentified number of equations"        if self.__numEquations==None: raise UndefinedPDEError("unidentified number of equations")
2409        if self.__numSolutions==None: raise UndefinedPDEError,"unidentified number of solutions"        if self.__numSolutions==None: raise UndefinedPDEError("unidentified number of solutions")
2410        # now we check the shape of the coefficient if numEquations and numSolutions are set:        # now we check the shape of the coefficient if numEquations and numSolutions are set:
2411        for i,d in coefficients.iteritems():        for i,d in list(coefficients.items()):
2412          try:          try:
2413             self.__COEFFICIENTS[i].setValue(self.getDomain(),             self.__COEFFICIENTS[i].setValue(self.getDomain(),
2414                       self.getNumEquations(),self.getNumSolutions(),                       self.getNumEquations(),self.getNumSolutions(),
2415                       self.reduceEquationOrder(),self.reduceSolutionOrder(),d)                       self.reduceEquationOrder(),self.reduceSolutionOrder(),d)
2416             self.alteredCoefficient(i)             self.alteredCoefficient(i)
2417          except IllegalCoefficientFunctionSpace,m:          except IllegalCoefficientFunctionSpace as m:
2418              # if the function space is wrong then we try the reduced version:              # if the function space is wrong then we try the reduced version:
2419              i_red=i+"_reduced"              i_red=i+"_reduced"
2420              if (not i_red in coefficients.keys()) and i_red in self.__COEFFICIENTS.keys():              if (not i_red in list(coefficients.keys())) and i_red in list(self.__COEFFICIENTS.keys()):
2421                  try:                  try:
2422                      self.__COEFFICIENTS[i_red].setValue(self.getDomain(),                      self.__COEFFICIENTS[i_red].setValue(self.getDomain(),
2423                                                        self.getNumEquations(),self.getNumSolutions(),                                                        self.getNumEquations(),self.getNumSolutions(),
2424                                                        self.reduceEquationOrder(),self.reduceSolutionOrder(),d)                                                        self.reduceEquationOrder(),self.reduceSolutionOrder(),d)
2425                      self.alteredCoefficient(i_red)                      self.alteredCoefficient(i_red)
2426                  except IllegalCoefficientValue,m:                  except IllegalCoefficientValue as m:
2427                      raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))                      raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
2428                  except IllegalCoefficientFunctionSpace,m:                  except IllegalCoefficientFunctionSpace as m:
2429                      raise IllegalCoefficientFunctionSpace("Coefficient %s:%s"%(i,m))                      raise IllegalCoefficientFunctionSpace("Coefficient %s:%s"%(i,m))
2430              else:              else:
2431                  raise IllegalCoefficientFunctionSpace("Coefficient %s:%s"%(i,m))                  raise IllegalCoefficientFunctionSpace("Coefficient %s:%s"%(i,m))
2432          except IllegalCoefficientValue,m:          except IllegalCoefficientValue as m:
2433             raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))             raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
2434        self.__altered_coefficients=True        self.__altered_coefficients=True
2435    
# Line 2740  class LinearPDE(LinearProblem): Line 2740  class LinearPDE(LinearProblem):
2740            mat,f=self.getSystem()            mat,f=self.getSystem()
2741            if self.isUsingLumping():            if self.isUsingLumping():
2742           if not util.inf(abs(mat)) > 0.:           if not util.inf(abs(mat)) > 0.:
2743           raise ZeroDivisionError,"Lumped mass matrix as zero entry (try order 1 elements or HRZ lumping)."           raise ZeroDivisionError("Lumped mass matrix as zero entry (try order 1 elements or HRZ lumping).")
2744               self.setSolution(f*1/mat)               self.setSolution(f*1/mat)
2745            else:            else:
2746               self.trace("PDE is resolved.")               self.trace("PDE is resolved.")
# Line 2760  class LinearPDE(LinearProblem): Line 2760  class LinearPDE(LinearProblem):
2760            if self.isUsingLumping():            if self.isUsingLumping():
2761                if not self.isOperatorValid():                if not self.isOperatorValid():
2762                   if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():                   if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():
2763                        raise TypeError,"Lumped matrix requires same order for equations and unknowns"                        raise TypeError("Lumped matrix requires same order for equations and unknowns")
2764                   if not self.getCoefficient("A").isEmpty():                   if not self.getCoefficient("A").isEmpty():
2765                        raise ValueError,"coefficient A in lumped matrix may not be present."                        raise ValueError("coefficient A in lumped matrix may not be present.")
2766                   if not self.getCoefficient("B").isEmpty():                   if not self.getCoefficient("B").isEmpty():
2767                        raise ValueError,"coefficient B in lumped matrix may not be present."                        raise ValueError("coefficient B in lumped matrix may not be present.")
2768                   if not self.getCoefficient("C").isEmpty():                   if not self.getCoefficient("C").isEmpty():
2769                        raise ValueError,"coefficient C in lumped matrix may not be present."                        raise ValueError("coefficient C in lumped matrix may not be present.")
2770                   if not self.getCoefficient("d_contact").isEmpty():                   if not self.getCoefficient("d_contact").isEmpty():
2771                        raise ValueError,"coefficient d_contact in lumped matrix may not be present."                        raise ValueError("coefficient d_contact in lumped matrix may not be present.")
2772                   if not self.getCoefficient("A_reduced").isEmpty():                   if not self.getCoefficient("A_reduced").isEmpty():
2773                        raise ValueError,"coefficient A_reduced in lumped matrix may not be present."                        raise ValueError("coefficient A_reduced in lumped matrix may not be present.")
2774                   if not self.getCoefficient("B_reduced").isEmpty():                   if not self.getCoefficient("B_reduced").isEmpty():
2775                        raise ValueError,"coefficient B_reduced in lumped matrix may not be present."                        raise ValueError("coefficient B_reduced in lumped matrix may not be present.")
2776                   if not self.getCoefficient("C_reduced").isEmpty():                   if not self.getCoefficient("C_reduced").isEmpty():
2777                        raise ValueError,"coefficient C_reduced in lumped matrix may not be present."                        raise ValueError("coefficient C_reduced in lumped matrix may not be present.")
2778                   if not self.getCoefficient("d_contact_reduced").isEmpty():                   if not self.getCoefficient("d_contact_reduced").isEmpty():
2779                        raise ValueError,"coefficient d_contact_reduced in lumped matrix may not be present."                        raise ValueError("coefficient d_contact_reduced in lumped matrix may not be present.")
2780                   D=self.getCoefficient("D")                   D=self.getCoefficient("D")
2781                   d=self.getCoefficient("d")                   d=self.getCoefficient("d")
2782                   D_reduced=self.getCoefficient("D_reduced")                   D_reduced=self.getCoefficient("D_reduced")
# Line 3835  class TransportPDE(LinearProblem): Line 3835  class TransportPDE(LinearProblem):
3835         :type value: large positive ``float``         :type value: large positive ``float``
3836         """         """
3837         if not value>0:         if not value>0:
3838           raise ValueError,"weighting factor needs to be positive."           raise ValueError("weighting factor needs to be positive.")
3839         self.__constraint_factor=value         self.__constraint_factor=value
3840         self.trace("Weighting factor for constraints is set to %e."%value)         self.trace("Weighting factor for constraints is set to %e."%value)
3841    
# Line 3863  class TransportPDE(LinearProblem): Line 3863  class TransportPDE(LinearProblem):
3863         if not dt == None:         if not dt == None:
3864        option_class=self.getSolverOptions()        option_class=self.getSolverOptions()
3865        if dt<=0:        if dt<=0:
3866            raise ValueError,"step size needs to be positive."            raise ValueError("step size needs to be positive.")
3867        if u0 == None:        if u0 == None:
3868            u0=self.getCurrentSolution()            u0=self.getCurrentSolution()
3869        else:        else:
3870            u0=util.interpolate(u0,self.getFunctionSpaceForSolution())            u0=util.interpolate(u0,self.getFunctionSpaceForSolution())
3871            if self.getNumSolutions() == 1:            if self.getNumSolutions() == 1:
3872          if u0.getShape()!=():          if u0.getShape()!=():
3873              raise ValueError,"Illegal shape %s of initial solution."%(u0.getShape(),)              raise ValueError("Illegal shape %s of initial solution."%(u0.getShape(),))
3874            else:            else:
3875              if u0.getShape()!=(self.getNumSolutions(),):              if u0.getShape()!=(self.getNumSolutions(),):
3876                raise ValueError,"Illegal shape %s of initial solution."%(u0.getShape(),)                raise ValueError("Illegal shape %s of initial solution."%(u0.getShape(),))
3877        self.setSolution(self.getOperator().solve(u0, self.getRightHandSide(),dt,option_class))        self.setSolution(self.getOperator().solve(u0, self.getRightHandSide(),dt,option_class))
3878        self.validSolution()        self.validSolution()
3879         return self.getCurrentSolution()         return self.getCurrentSolution()
# Line 3889  class TransportPDE(LinearProblem): Line 3889  class TransportPDE(LinearProblem):
3889         u2=util.interpolate(u,self.getFunctionSpaceForSolution())         u2=util.interpolate(u,self.getFunctionSpaceForSolution())
3890         if self.getNumSolutions() == 1:         if self.getNumSolutions() == 1:
3891            if u2.getShape()!=():            if u2.getShape()!=():
3892                raise ValueError,"Illegal shape %s of initial solution."%(u2.getShape(),)                raise ValueError("Illegal shape %s of initial solution."%(u2.getShape(),))
3893         else:         else:
3894            if u2.getShape()!=(self.getNumSolutions(),):            if u2.getShape()!=(self.getNumSolutions(),):
3895                raise ValueError,"Illegal shape %s of initial solution."%(u2.getShape(),)                raise ValueError("Illegal shape %s of initial solution."%(u2.getShape(),))
3896         self.setSolution(u2,validate=False)         self.setSolution(u2,validate=False)
3897    
3898    

Legend:
Removed from v.3788  
changed lines
  Added in v.3789

  ViewVC Help
Powered by ViewVC 1.1.26