/[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 3851 by gross, Mon Feb 13 06:22:06 2012 UTC revision 3852 by jfenwick, Thu Mar 1 05:34:52 2012 UTC
# Line 202  class SolverOptions(object): Line 202  class SolverOptions(object):
202          self.setUsePanel()          self.setUsePanel()
203          self.setDiagonalDominanceThreshold()          self.setDiagonalDominanceThreshold()
204          self.setAMGInterpolation()          self.setAMGInterpolation()
205      self.setCycleType()          self.setCycleType()
206      self.setODESolver()          self.setODESolver()
207                    
208    
209      def __str__(self):      def __str__(self):
# Line 252  class SolverOptions(object): Line 252  class SolverOptions(object):
252                  out+="\nCoarsening threshold = %e"%self.getMinCoarseMatrixSize()                  out+="\nCoarsening threshold = %e"%self.getMinCoarseMatrixSize()
253                  out+="\nMinimum size of the coarsest level matrix = %e"%self.getCoarseningThreshold()                  out+="\nMinimum size of the coarsest level matrix = %e"%self.getCoarseningThreshold()
254                  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())
255          if self.getPreconditioner() == self.BOOMERAMG:              if self.getPreconditioner() == self.BOOMERAMG:
256          out+="\nMaximum number of levels = %s"%self.getLevelMax()                  out+="\nMaximum number of levels = %s"%self.getLevelMax()
257          out+="\nCoarsening threshold = %e"%self.getCoarseningThreshold()                  out+="\nCoarsening threshold = %e"%self.getCoarseningThreshold()
258          out+="\nThreshold for diagonal dominant rows = %s"%(setDiagonalDominanceThreshold(),)                  out+="\nThreshold for diagonal dominant rows = %s"%(setDiagonalDominanceThreshold(),)
259          out+="\nCoarsening method = %s"%self.getName(self.getCoarsening())                  out+="\nCoarsening method = %s"%self.getName(self.getCoarsening())
260          out+="\nV-cycle (1) or W-cyle (2) = %s"%self.getCycleType()                  out+="\nV-cycle (1) or W-cyle (2) = %s"%self.getCycleType()
261          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())
262          out+="\nSmoother = %s"%self.getName(self.getSmoother())                  out+="\nSmoother = %s"%self.getName(self.getSmoother())
263          if self.getPreconditioner() == self.GAUSS_SEIDEL:              if self.getPreconditioner() == self.GAUSS_SEIDEL:
264                  out+="\nNumber of sweeps = %s"%self.getNumSweeps()                  out+="\nNumber of sweeps = %s"%self.getNumSweeps()
265              if self.getPreconditioner() == self.ILUT:              if self.getPreconditioner() == self.ILUT:
266                  out+="\nDrop tolerance = %e"%self.getDropTolerance()                  out+="\nDrop tolerance = %e"%self.getDropTolerance()
# Line 317  class SolverOptions(object): Line 317  class SolverOptions(object):
317          if key == self.CLASSIC_INTERPOLATION_WITH_FF_COUPLING: return "CLASSIC_INTERPOLATION_WITH_FF"          if key == self.CLASSIC_INTERPOLATION_WITH_FF_COUPLING: return "CLASSIC_INTERPOLATION_WITH_FF"
318          if key == self.CLASSIC_INTERPOLATION: return "CLASSIC_INTERPOLATION"          if key == self.CLASSIC_INTERPOLATION: return "CLASSIC_INTERPOLATION"
319          if key == self.DIRECT_INTERPOLATION: return "DIRECT_INTERPOLATION"          if key == self.DIRECT_INTERPOLATION: return "DIRECT_INTERPOLATION"
320      if key == self.BOOMERAMG: return "BOOMERAMG"          if key == self.BOOMERAMG: return "BOOMERAMG"
321      if key == self.CIJP_FIXED_RANDOM_COARSENING: return "CIJP_FIXED_RANDOM_COARSENING"          if key == self.CIJP_FIXED_RANDOM_COARSENING: return "CIJP_FIXED_RANDOM_COARSENING"
322      if key == self.CIJP_COARSENING: return "CIJP_COARSENING"          if key == self.CIJP_COARSENING: return "CIJP_COARSENING"
323      if key == self.FALGOUT_COARSENING: return "FALGOUT_COARSENING"          if key == self.FALGOUT_COARSENING: return "FALGOUT_COARSENING"
324      if key == self.PMIS_COARSENING: return "PMIS_COARSENING"          if key == self.PMIS_COARSENING: return "PMIS_COARSENING"
325      if key == self.HMIS_COARSENING: return "HMIS_COARSENING"          if key == self.HMIS_COARSENING: return "HMIS_COARSENING"
326      if key == self.LINEAR_CRANK_NICOLSON: return "LINEAR_CRANK_NICOLSON"          if key == self.LINEAR_CRANK_NICOLSON: return "LINEAR_CRANK_NICOLSON"
327          if key == self.CRANK_NICOLSON: return "CRANK_NICOLSON"          if key == self.CRANK_NICOLSON: return "CRANK_NICOLSON"
328          if key == self.BACKWARD_EULER: return "BACKWARD_EULER"          if key == self.BACKWARD_EULER: return "BACKWARD_EULER"
329    
# Line 386  class SolverOptions(object): Line 386  class SolverOptions(object):
386          if name == "time_step_backtracking_used":          if name == "time_step_backtracking_used":
387              self.__time_step_backtracking_used = (value == True)              self.__time_step_backtracking_used = (value == True)
388          if name == "coarse_level_sparsity":          if name == "coarse_level_sparsity":
389             self.__coarse_level_sparsity = value              self.__coarse_level_sparsity = value
390      if name == "num_coarse_unknowns":          if name == "num_coarse_unknowns":
391             self.__num_coarse_unknowns= value              self.__num_coarse_unknowns= value
392      def getDiagnostics(self, name):      def getDiagnostics(self, name):
393          """          """
394          Returns the diagnostic information ``name``. Possible values are:          Returns the diagnostic information ``name``. Possible values are:
# Line 448  class SolverOptions(object): Line 448  class SolverOptions(object):
448          :param method: selects the coarsening method .          :param method: selects the coarsening method .
449          :type method: in {SolverOptions.DEFAULT}, `SolverOptions.YAIR_SHAPIRA_COARSENING`,  `SolverOptions.RUGE_STUEBEN_COARSENING`, `SolverOptions.AGGREGATION_COARSENING`, `SolverOptions.CIJP_FIXED_RANDOM_COARSENING`, `SolverOptions.CIJP_COARSENING`, `SolverOptions.FALGOUT_COARSENING`, `SolverOptions.PMIS_COARSENING`, `SolverOptions.HMIS_COARSENING`          :type method: in {SolverOptions.DEFAULT}, `SolverOptions.YAIR_SHAPIRA_COARSENING`,  `SolverOptions.RUGE_STUEBEN_COARSENING`, `SolverOptions.AGGREGATION_COARSENING`, `SolverOptions.CIJP_FIXED_RANDOM_COARSENING`, `SolverOptions.CIJP_COARSENING`, `SolverOptions.FALGOUT_COARSENING`, `SolverOptions.PMIS_COARSENING`, `SolverOptions.HMIS_COARSENING`
450          """          """
451      if method==None: method=0          if method==None: method=0
452          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,]:
453               raise ValueError("unknown coarsening method %s"%method)               raise ValueError("unknown coarsening method %s"%method)
454          self.__coarsening=method          self.__coarsening=method
# Line 468  class SolverOptions(object): Line 468  class SolverOptions(object):
468          :param size: minumum size of the coarsest level matrix .          :param size: minumum size of the coarsest level matrix .
469          :type size: positive ``int`` or ``None``          :type size: positive ``int`` or ``None``
470          """          """
471      if size==None: size=500          if size==None: size=500
472          size=int(size)          size=int(size)
473          if size<0:          if size<0:
474             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.")
# Line 492  class SolverOptions(object): Line 492  class SolverOptions(object):
492                                      `SolverOptions.NO_PRECONDITIONER`                                      `SolverOptions.NO_PRECONDITIONER`
493          :note: Not all packages support all preconditioner. It can be assumed that a package makes a reasonable choice if it encounters an unknown preconditioner.          :note: Not all packages support all preconditioner. It can be assumed that a package makes a reasonable choice if it encounters an unknown preconditioner.
494          """          """
495      if preconditioner==None: preconditioner=10          if preconditioner==None: preconditioner=10
496          if not preconditioner in [  SolverOptions.ILU0, SolverOptions.ILUT, SolverOptions.JACOBI,          if not preconditioner in [  SolverOptions.ILU0, SolverOptions.ILUT, SolverOptions.JACOBI,
497                                      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,
498                                      SolverOptions.NO_PRECONDITIONER] :                                      SolverOptions.NO_PRECONDITIONER] :
# Line 515  class SolverOptions(object): Line 515  class SolverOptions(object):
515          :type smoother: in `SolverOptions.JACOBI`, `SolverOptions.GAUSS_SEIDEL`          :type smoother: in `SolverOptions.JACOBI`, `SolverOptions.GAUSS_SEIDEL`
516          :note: Not all packages support all smoothers. It can be assumed that a package makes a reasonable choice if it encounters an unknown smoother.          :note: Not all packages support all smoothers. It can be assumed that a package makes a reasonable choice if it encounters an unknown smoother.
517          """          """
518      if smoother==None: smoother=28          if smoother==None: smoother=28
519          if not smoother in [ SolverOptions.JACOBI, SolverOptions.GAUSS_SEIDEL ] :          if not smoother in [ SolverOptions.JACOBI, SolverOptions.GAUSS_SEIDEL ] :
520               raise ValueError("unknown smoother %s"%smoother)               raise ValueError("unknown smoother %s"%smoother)
521          self.__smoother=smoother              self.__smoother=smoother    
# Line 540  class SolverOptions(object): Line 540  class SolverOptions(object):
540                          `SolverOptions.NONLINEAR_GMRES`, `SolverOptions.TFQMR`, `SolverOptions.MINRES`                          `SolverOptions.NONLINEAR_GMRES`, `SolverOptions.TFQMR`, `SolverOptions.MINRES`
541          :note: Not all packages support all solvers. It can be assumed that a package makes a reasonable choice if it encounters an unknown solver method.          :note: Not all packages support all solvers. It can be assumed that a package makes a reasonable choice if it encounters an unknown solver method.
542          """          """
543      if method==None: method=0          if method==None: method=0
544          if not method in [ SolverOptions.DEFAULT, SolverOptions.DIRECT, SolverOptions.CHOLEVSKY, SolverOptions.PCG,          if not method in [ SolverOptions.DEFAULT, SolverOptions.DIRECT, SolverOptions.CHOLEVSKY, SolverOptions.PCG,
545                             SolverOptions.CR, SolverOptions.CGS, SolverOptions.BICGSTAB,                             SolverOptions.CR, SolverOptions.CGS, SolverOptions.BICGSTAB,
546                             SolverOptions.GMRES, SolverOptions.PRES20, SolverOptions.ROWSUM_LUMPING, SolverOptions.HRZ_LUMPING,                             SolverOptions.GMRES, SolverOptions.PRES20, SolverOptions.ROWSUM_LUMPING, SolverOptions.HRZ_LUMPING,
# Line 567  class SolverOptions(object): Line 567  class SolverOptions(object):
567          :type package: in `SolverOptions.DEFAULT`, `SolverOptions.PASO`, `SolverOptions.SUPER_LU`, `SolverOptions.PASTIX`, `SolverOptions.MKL`, `SolverOptions.UMFPACK`, `SolverOptions.TRILINOS`          :type package: in `SolverOptions.DEFAULT`, `SolverOptions.PASO`, `SolverOptions.SUPER_LU`, `SolverOptions.PASTIX`, `SolverOptions.MKL`, `SolverOptions.UMFPACK`, `SolverOptions.TRILINOS`
568          :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.
569          """          """
570      if package==None: package=0          if package==None: package=0
571          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]:
572               raise ValueError("unknown solver package %s"%package)               raise ValueError("unknown solver package %s"%package)
573          self.__package=package          self.__package=package
# Line 623  class SolverOptions(object): Line 623  class SolverOptions(object):
623          else:          else:
624              return self.__restart              return self.__restart
625      def _getRestartForC(self):      def _getRestartForC(self):
626          r=self.getRestart()              r=self.getRestart()
627          if r == None:              if r == None:
628              return -1                 return -1
629              else:              else:
630              return r                 return r
631        
632      def setDiagonalDominanceThreshold(self,value=0.5):      def setDiagonalDominanceThreshold(self,value=0.5):
633          """          """
# Line 638  class SolverOptions(object): Line 638  class SolverOptions(object):
638          """          """
639          value=float(value)          value=float(value)
640          if value<0 or value>1.:          if value<0 or value>1.:
641         raise ValueError("Diagonal dominance threshold must be between 0 and 1.")             raise ValueError("Diagonal dominance threshold must be between 0 and 1.")
642      self.__diagonal_dominance_threshold=value          self.__diagonal_dominance_threshold=value
643                    
644      def getDiagonalDominanceThreshold(self):      def getDiagonalDominanceThreshold(self):
645          """          """
# Line 1097  class SolverOptions(object): Line 1097  class SolverOptions(object):
1097         """         """
1098         sparsity=float(sparsity)         sparsity=float(sparsity)
1099         if sparsity<0:         if sparsity<0:
1100       raise ValueError("sparsity must be non-negative.")            raise ValueError("sparsity must be non-negative.")
1101         if sparsity>1:         if sparsity>1:
1102            raise ValueError("sparsity must be less than 1.")            raise ValueError("sparsity must be less than 1.")
1103         self.__min_sparsity=sparsity         self.__min_sparsity=sparsity
# Line 1122  class SolverOptions(object): Line 1122  class SolverOptions(object):
1122        """        """
1123        refinements=int(refinements)        refinements=int(refinements)
1124        if refinements<0:        if refinements<0:
1125       raise ValueError("number of refinements must be non-negative.")          raise ValueError("number of refinements must be non-negative.")
1126        self.__refinements=refinements        self.__refinements=refinements
1127        
1128      def getNumRefinements(self):      def getNumRefinements(self):
# Line 1142  class SolverOptions(object): Line 1142  class SolverOptions(object):
1142        """        """
1143        refinements=int(refinements)        refinements=int(refinements)
1144        if refinements<0:        if refinements<0:
1145       raise ValueError("number of refinements must be non-negative.")           raise ValueError("number of refinements must be non-negative.")
1146        self.__coarse_refinements=refinements        self.__coarse_refinements=refinements
1147        
1148      def getNumCoarseMatrixRefinements(self):      def getNumCoarseMatrixRefinements(self):
# Line 1194  class SolverOptions(object): Line 1194  class SolverOptions(object):
1194          :param method: key of the interpolation method to be used.          :param method: key of the interpolation method to be used.
1195          :type method: in `SolverOptions.CLASSIC_INTERPOLATION_WITH_FF_COUPLING`, `SolverOptions.CLASSIC_INTERPOLATION', `SolverOptions.DIRECT_INTERPOLATION`          :type method: in `SolverOptions.CLASSIC_INTERPOLATION_WITH_FF_COUPLING`, `SolverOptions.CLASSIC_INTERPOLATION', `SolverOptions.DIRECT_INTERPOLATION`
1196          """          """
1197      if method==None: method=self.DIRECT_INTERPOLATION          if method==None: method=self.DIRECT_INTERPOLATION
1198          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 ]:
1199               raise ValueError("unknown AMG interpolation method %s"%method)               raise ValueError("unknown AMG interpolation method %s"%method)
1200          self.__amg_interpolation_method=method          self.__amg_interpolation_method=method
# Line 1214  class SolverOptions(object): Line 1214  class SolverOptions(object):
1214          :param method: key of the ODE solver method to be used.          :param method: key of the ODE solver method to be used.
1215          :type method: in `SolverOptions.CRANK_NICOLSON, `SolverOptions.BACKWARD_EULER', `SolverOptions.LINEAR_CRANK_NICOLSON          :type method: in `SolverOptions.CRANK_NICOLSON, `SolverOptions.BACKWARD_EULER', `SolverOptions.LINEAR_CRANK_NICOLSON
1216          """          """
1217      if method==None: method=self.LINEAR_CRANK_NICOLSON          if method==None: method=self.LINEAR_CRANK_NICOLSON
1218          if not method in [ SolverOptions.CRANK_NICOLSON, SolverOptions.BACKWARD_EULER, SolverOptions.LINEAR_CRANK_NICOLSON ]:          if not method in [ SolverOptions.CRANK_NICOLSON, SolverOptions.BACKWARD_EULER, SolverOptions.LINEAR_CRANK_NICOLSON ]:
1219               raise ValueError("unknown ODE solver method %s"%method)               raise ValueError("unknown ODE solver method %s"%method)
1220          self.__ode_solver=method          self.__ode_solver=method
# Line 2232  class LinearProblem(object): Line 2232  class LinearProblem(object):
2232         if not self.getDomainStatus()==self.getSystemStatus(): self.invalidateSolution()         if not self.getDomainStatus()==self.getSystemStatus(): self.invalidateSolution()
2233         if self.__solution_rtol>self.getSolverOptions().getTolerance() or \         if self.__solution_rtol>self.getSolverOptions().getTolerance() or \
2234            self.__solution_atol>self.getSolverOptions().getAbsoluteTolerance():            self.__solution_atol>self.getSolverOptions().getAbsoluteTolerance():
2235           self.invalidateSolution()                self.invalidateSolution()  
2236         return self.__is_solution_valid         return self.__is_solution_valid
2237    
2238     def validOperator(self):     def validOperator(self):
# Line 2357  class LinearProblem(object): Line 2357  class LinearProblem(object):
2357         defined by the solver options         defined by the solver options
2358         """         """
2359         if validate:         if validate:
2360        self.__solution_rtol=self.getSolverOptions().getTolerance()            self.__solution_rtol=self.getSolverOptions().getTolerance()
2361        self.__solution_atol=self.getSolverOptions().getAbsoluteTolerance()            self.__solution_atol=self.getSolverOptions().getAbsoluteTolerance()
2362        self.validSolution()            self.validSolution()
2363         self.__solution=u         self.__solution=u
2364     def getCurrentSolution(self):     def getCurrentSolution(self):
2365         """         """
# Line 2395  class LinearProblem(object): Line 2395  class LinearProblem(object):
2395                 self.__operator=self.createSolution()                 self.__operator=self.createSolution()
2396             else:             else:
2397                 self.__operator=self.createOperator()                 self.__operator=self.createOperator()
2398         self.__operator_type=self.getRequiredOperatorType()             self.__operator_type=self.getRequiredOperatorType()
2399         else:         else:
2400             if self.isUsingLumping():             if self.isUsingLumping():
2401                 self.__operator.setToZero()                 self.__operator.setToZero()
# Line 2404  class LinearProblem(object): Line 2404  class LinearProblem(object):
2404                     self.__operator.resetValues()                     self.__operator.resetValues()
2405                 else:                 else:
2406                     self.__operator=self.createOperator()                     self.__operator=self.createOperator()
2407                 self.__operator_type=self.getRequiredOperatorType()                     self.__operator_type=self.getRequiredOperatorType()
2408             self.trace("Operator reset to zero")             self.trace("Operator reset to zero")
2409    
2410     def getCurrentOperator(self):     def getCurrentOperator(self):
# Line 2719  class LinearPDE(LinearProblem): Line 2719  class LinearPDE(LinearProblem):
2719        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.
2720        """        """
2721        if self.isUsingLumping():        if self.isUsingLumping():
2722       return "__ESCRIPT_DATA"           return "__ESCRIPT_DATA"
2723        else:        else:
2724       solver_options=self.getSolverOptions()           solver_options=self.getSolverOptions()
2725       return self.getDomain().getSystemMatrixTypeId(solver_options.getSolverMethod(), solver_options.getPreconditioner(),solver_options.getPackage(), solver_options.isSymmetric())           return self.getDomain().getSystemMatrixTypeId(solver_options.getSolverMethod(), solver_options.getPreconditioner(),solver_options.getPackage(), solver_options.isSymmetric())
2726    
2727     def checkSymmetry(self,verbose=True):     def checkSymmetry(self,verbose=True):
2728        """        """
# Line 2774  class LinearPDE(LinearProblem): Line 2774  class LinearPDE(LinearProblem):
2774         if not self.isSolutionValid():         if not self.isSolutionValid():
2775            mat,f=self.getSystem()            mat,f=self.getSystem()
2776            if self.isUsingLumping():            if self.isUsingLumping():
2777           if not util.inf(abs(mat)) > 0.:               if not util.inf(abs(mat)) > 0.:
2778           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).")
2779               self.setSolution(f*1/mat)               self.setSolution(f*1/mat)
2780            else:            else:
2781               self.trace("PDE is resolved.")               self.trace("PDE is resolved.")
# Line 2860  class LinearPDE(LinearProblem): Line 2860  class LinearPDE(LinearProblem):
2860                   self.resetOperator()                   self.resetOperator()
2861                   operator=self.getCurrentOperator()                   operator=self.getCurrentOperator()
2862                   if hasattr(self.getDomain(), "addPDEToLumpedSystem") :                   if hasattr(self.getDomain(), "addPDEToLumpedSystem") :
2863              hrz_lumping=( self.getSolverOptions().getSolverMethod() ==  SolverOptions.HRZ_LUMPING )                      hrz_lumping=( self.getSolverOptions().getSolverMethod() ==  SolverOptions.HRZ_LUMPING )
2864              self.getDomain().addPDEToLumpedSystem(operator, D_times_e, d_times_e, d_dirac_times_e,  hrz_lumping )                      self.getDomain().addPDEToLumpedSystem(operator, D_times_e, d_times_e, d_dirac_times_e,  hrz_lumping )
2865              self.getDomain().addPDEToLumpedSystem(operator, D_reduced_times_e, d_reduced_times_e, escript.Data(), hrz_lumping)                      self.getDomain().addPDEToLumpedSystem(operator, D_reduced_times_e, d_reduced_times_e, escript.Data(), hrz_lumping)
2866                   else:                   else:
2867                      self.getDomain().addPDEToRHS(operator, \                      self.getDomain().addPDEToRHS(operator, \
2868                                                   escript.Data(), \                                                   escript.Data(), \
# Line 3864  class TransportPDE(LinearProblem): Line 3864  class TransportPDE(LinearProblem):
3864         :rtype: `Data`         :rtype: `Data`
3865         """         """
3866         if not dt == None:         if not dt == None:
3867        option_class=self.getSolverOptions()            option_class=self.getSolverOptions()
3868        if dt<=0:            if dt<=0:
3869            raise ValueError("step size needs to be positive.")                raise ValueError("step size needs to be positive.")
3870        if u0 == None:            if u0 == None:
3871            u0=self.getCurrentSolution()                u0=self.getCurrentSolution()
3872        else:            else:
3873            u0=util.interpolate(u0,self.getFunctionSpaceForSolution())                u0=util.interpolate(u0,self.getFunctionSpaceForSolution())
3874            if self.getNumSolutions() == 1:                if self.getNumSolutions() == 1:
3875          if u0.getShape()!=():                  if u0.getShape()!=():
3876              raise ValueError("Illegal shape %s of initial solution."%(u0.getShape(),))                    raise ValueError("Illegal shape %s of initial solution."%(u0.getShape(),))
3877            else:                else:
3878              if u0.getShape()!=(self.getNumSolutions(),):                   if u0.getShape()!=(self.getNumSolutions(),):
3879                raise ValueError("Illegal shape %s of initial solution."%(u0.getShape(),))                     raise ValueError("Illegal shape %s of initial solution."%(u0.getShape(),))
3880        self.setSolution(self.getOperator().solve(u0, self.getRightHandSide(),dt,option_class))            self.setSolution(self.getOperator().solve(u0, self.getRightHandSide(),dt,option_class))
3881        self.validSolution()            self.validSolution()
3882         return self.getCurrentSolution()         return self.getCurrentSolution()
3883    
3884     def setInitialSolution(self,u):     def setInitialSolution(self,u):

Legend:
Removed from v.3851  
changed lines
  Added in v.3852

  ViewVC Help
Powered by ViewVC 1.1.26