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

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

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

revision 1557 by artak, Mon May 19 04:44:27 2008 UTC revision 1737 by gross, Fri Aug 29 03:01:29 2008 UTC
# Line 474  class IterationHistory(object): Line 474  class IterationHistory(object):
474    
475         """         """
476         self.history.append(norm_r)         self.history.append(norm_r)
477         if self.verbose: print "iter: #s:  inner(rhat,r) = #e"#(len(self.history)-1, self.history[-1])         if self.verbose: print "iter: %s:  inner(rhat,r) = %e"#(len(self.history)-1, self.history[-1])
478         return self.history[-1]<=self.tolerance * self.history[0]         return self.history[-1]<=self.tolerance * self.history[0]
479    
480     def stoppingcriterium2(self,norm_r,norm_b,solver="GMRES",TOL=None):     def stoppingcriterium2(self,norm_r,norm_b,solver="GMRES",TOL=None):
# Line 493  class IterationHistory(object): Line 493  class IterationHistory(object):
493         if TOL==None:         if TOL==None:
494            TOL=self.tolerance            TOL=self.tolerance
495         self.history.append(norm_r)         self.history.append(norm_r)
496         if self.verbose: print "iter: #s:  norm(r) = #e"#(len(self.history)-1, self.history[-1])         if self.verbose: print "iter: %s:  norm(r) = %e"#(len(self.history)-1, self.history[-1])
497         return self.history[-1]<=TOL * norm_b         return self.history[-1]<=TOL * norm_b
498    
499  def PCG(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):  def PCG(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
# Line 1044  def NewtonGMRES(b, Aprod, Msolve, biline Line 1044  def NewtonGMRES(b, Aprod, Msolve, biline
1044              etanew=gamma*rat*rat              etanew=gamma*rat*rat
1045              if gamma*etaold*etaold > .1 :              if gamma*etaold*etaold > .1 :
1046                  etanew=max(etanew,gamma*etaold*etaold)                  etanew=max(etanew,gamma*etaold*etaold)
               
1047              etamax=min(etanew,etamax)              etamax=min(etanew,etamax)
1048              etamax=max(etamax,.5*stop_tol/fnrm)              etamax=max(etamax,.5*stop_tol/fnrm)
   
1049      return x      return x
1050    
1051  def TFQMR(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):  def TFQMR(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
# Line 1514  class HomogeneousSaddlePointProblem(obje Line 1512  class HomogeneousSaddlePointProblem(obje
1512                  if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter                  if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter
1513                  p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p)                  p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p)
1514              u=r[0]                u=r[0]  
1515                  print "DIFF=",util.integrate(p)                  # print "DIFF=",util.integrate(p)
1516    
1517                print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)                # print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)
1518    
1519            return u,p            return u,p
1520    
# Line 1592  class SaddlePointProblem(object): Line 1590  class SaddlePointProblem(object):
1590         @type verbose: C{bool}         @type verbose: C{bool}
1591         @note: this method may be overwritten by a particular saddle point problem         @note: this method may be overwritten by a particular saddle point problem
1592         """         """
1593           print "SaddlePointProblem should not be used anymore!"
1594         if not isinstance(verbose,bool):         if not isinstance(verbose,bool):
1595              raise TypeError("verbose needs to be of type bool.")              raise TypeError("verbose needs to be of type bool.")
1596         self.__verbose=verbose         self.__verbose=verbose
# Line 1604  class SaddlePointProblem(object): Line 1603  class SaddlePointProblem(object):
1603         @param text: a text message         @param text: a text message
1604         @type text: C{str}         @type text: C{str}
1605         """         """
1606         if self.__verbose: print "#s: #s"%(str(self),text)         if self.__verbose: print "%s: %s"%(str(self),text)
1607    
1608     def solve_f(self,u,p,tol=1.e-8):     def solve_f(self,u,p,tol=1.e-8):
1609         """         """
# Line 1796  class SaddlePointProblem(object): Line 1795  class SaddlePointProblem(object):
1795  #  #
1796  #      return u,p  #      return u,p
1797                        
1798  def MaskFromBoundaryTag(function_space,*tags):  def MaskFromBoundaryTag(domain,*tags):
1799     """     """
1800     create a mask on the given function space which one for samples     create a mask on the Solution(domain) function space which one for samples
1801     that touch the boundary tagged by tags.     that touch the boundary tagged by tags.
1802    
1803     usage: m=MaskFromBoundaryTag(Solution(domain),"left", "right")     usage: m=MaskFromBoundaryTag(domain,"left", "right")
1804    
1805     @param function_space: a given function space     @param domain: a given domain
1806     @type function_space: L{escript.FunctionSpace}     @type domain: L{escript.Domain}
1807     @param tags: boundray tags     @param tags: boundray tags
1808     @type tags: C{str}     @type tags: C{str}
1809     @return: a mask which marks samples used by C{function_space} that are touching the     @return: a mask which marks samples that are touching the boundary tagged by any of the given tags.
             boundary tagged by any of the given tags.  
1810     @rtype: L{escript.Data} of rank 0     @rtype: L{escript.Data} of rank 0
1811     """     """
1812     pde=linearPDEs.LinearPDE(function_space.getDomain(),numEquations=1, numSolutions=1)     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1813     d=escript.Scalar(0.,escript.FunctionOnBoundary(function_space.getDomain()))     d=escript.Scalar(0.,escript.FunctionOnBoundary(function_space.getDomain()))
1814     for t in tags: d.setTaggedValue(t,1.)     for t in tags: d.setTaggedValue(t,1.)
1815     pde.setValue(y=d)     pde.setValue(y=d)
1816     out=util.whereNonZero(pde.getRightHandSide())     return util.whereNonZero(pde.getRightHandSide())
    if out.getFunctionSpace() == function_space:  
       return out  
    else:  
       return util.whereNonZero(util.interpolate(out,function_space))  
   
   
   

Legend:
Removed from v.1557  
changed lines
  Added in v.1737

  ViewVC Help
Powered by ViewVC 1.1.26