/[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 3676 by jfenwick, Thu Nov 17 03:52:25 2011 UTC revision 3771 by jfenwick, Wed Jan 18 02:30:48 2012 UTC
# Line 39  Currently includes: Line 39  Currently includes:
39  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
40    
41    
42  import escript  from . import escript
43  import linearPDEs  from . import linearPDEs
44  import numpy  import numpy
45  import util  from . import util
46  import math  import math
47    
48  class TimeIntegrationManager:  class TimeIntegrationManager:
# Line 66  class TimeIntegrationManager: Line 66  class TimeIntegrationManager:
66       Sets up the value manager where ``inital_values`` are the initial values       Sets up the value manager where ``inital_values`` are the initial values
67       and p is the order used for extrapolation.       and p is the order used for extrapolation.
68       """       """
69       if kwargs.has_key("p"):       if "p" in kwargs:
70              self.__p=kwargs["p"]              self.__p=kwargs["p"]
71       else:       else:
72              self.__p=1              self.__p=1
73       if kwargs.has_key("time"):       if "time" in kwargs:
74              self.__t=kwargs["time"]              self.__t=kwargs["time"]
75       else:       else:
76              self.__t=0.              self.__t=0.
# Line 299  class NoPDE: Line 299  class NoPDE:
299           """           """
300           if self.__u==None:           if self.__u==None:
301              if self.__D==None:              if self.__D==None:
302                 raise ValueError,"coefficient D is undefined"                 raise ValueError("coefficient D is undefined")
303              D=escript.Data(self.__D,self.__function_space)              D=escript.Data(self.__D,self.__function_space)
304              if D.getRank()>1:              if D.getRank()>1:
305                 raise ValueError,"coefficient D must have rank 0 or 1"                 raise ValueError("coefficient D must have rank 0 or 1")
306              if self.__Y==None:              if self.__Y==None:
307                 self.__u=escript.Data(0.,D.getShape(),self.__function_space)                 self.__u=escript.Data(0.,D.getShape(),self.__function_space)
308              else:              else:
# Line 340  class Locator: Line 340  class Locator:
340         iterative=False         iterative=False
341         if isinstance(x, list):         if isinstance(x, list):
342             if len(x)==0:             if len(x)==0:
343                raise ValueError, "At least one point must be given."                raise ValueError("At least one point must be given.")
344             try:             try:
345               iter(x[0])               iter(x[0])
346               iterative=True               iterative=True
# Line 452  def getInfLocator(arg): Line 452  def getInfLocator(arg):
452      Return a Locator for a point with the inf value over all arg.      Return a Locator for a point with the inf value over all arg.
453      """      """
454      if not isinstance(arg, escript.Data):      if not isinstance(arg, escript.Data):
455      raise TypeError,"getInfLocator: Unknown argument type."      raise TypeError("getInfLocator: Unknown argument type.")
456      a_inf=util.inf(arg)      a_inf=util.inf(arg)
457      loc=util.length(arg-a_inf).minGlobalDataPoint() # This gives us the location but not coords      loc=util.length(arg-a_inf).minGlobalDataPoint() # This gives us the location but not coords
458      x=arg.getFunctionSpace().getX()      x=arg.getFunctionSpace().getX()
# Line 464  def getSupLocator(arg): Line 464  def getSupLocator(arg):
464      Return a Locator for a point with the sup value over all arg.      Return a Locator for a point with the sup value over all arg.
465      """      """
466      if not isinstance(arg, escript.Data):      if not isinstance(arg, escript.Data):
467      raise TypeError,"getInfLocator: Unknown argument type."      raise TypeError("getInfLocator: Unknown argument type.")
468      a_inf=util.sup(arg)      a_inf=util.sup(arg)
469      loc=util.length(arg-a_inf).minGlobalDataPoint() # This gives us the location but not coords      loc=util.length(arg-a_inf).minGlobalDataPoint() # This gives us the location but not coords
470      x=arg.getFunctionSpace().getX()      x=arg.getFunctionSpace().getX()
# Line 563  def PCG(r, Aprod, x, Msolve, bilinearfor Line 563  def PCG(r, Aprod, x, Msolve, bilinearfor
563     rhat=Msolve(r)     rhat=Msolve(r)
564     d = rhat     d = rhat
565     rhat_dot_r = bilinearform(rhat, r)     rhat_dot_r = bilinearform(rhat, r)
566     if rhat_dot_r<0: raise NegativeNorm,"negative norm."     if rhat_dot_r<0: raise NegativeNorm("negative norm.")
567     norm_r0=math.sqrt(rhat_dot_r)     norm_r0=math.sqrt(rhat_dot_r)
568     atol2=atol+rtol*norm_r0     atol2=atol+rtol*norm_r0
569     if atol2<=0:     if atol2<=0:
570        raise ValueError,"Non-positive tolarance."        raise ValueError("Non-positive tolarance.")
571     atol2=max(atol2, 100. * util.EPSILON * norm_r0)     atol2=max(atol2, 100. * util.EPSILON * norm_r0)
572    
573     if verbose: print "PCG: initial residual norm = %e (absolute tolerance = %e)"%(norm_r0, atol2)     if verbose: print(("PCG: initial residual norm = %e (absolute tolerance = %e)"%(norm_r0, atol2)))
574    
575    
576     while not math.sqrt(rhat_dot_r) <= atol2:     while not math.sqrt(rhat_dot_r) <= atol2:
577         iter+=1         iter+=1
578         if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max         if iter  >= iter_max: raise MaxIterReached("maximum number of %s steps reached."%iter_max)
579    
580         q=Aprod(d)         q=Aprod(d)
581         alpha = rhat_dot_r / bilinearform(d, q)         alpha = rhat_dot_r / bilinearform(d, q)
# Line 591  def PCG(r, Aprod, x, Msolve, bilinearfor Line 591  def PCG(r, Aprod, x, Msolve, bilinearfor
591         d=rhat         d=rhat
592    
593         rhat_dot_r = rhat_dot_r_new         rhat_dot_r = rhat_dot_r_new
594         if rhat_dot_r<0: raise NegativeNorm,"negative norm."         if rhat_dot_r<0: raise NegativeNorm("negative norm.")
595         if verbose: print "PCG: iteration step %s: residual norm = %e"%(iter, math.sqrt(rhat_dot_r))         if verbose: print(("PCG: iteration step %s: residual norm = %e"%(iter, math.sqrt(rhat_dot_r))))
596     if verbose: print "PCG: tolerance reached after %s steps."%iter     if verbose: print(("PCG: tolerance reached after %s steps."%iter))
597     return x,r,math.sqrt(rhat_dot_r)     return x,r,math.sqrt(rhat_dot_r)
598    
599  class Defect(object):  class Defect(object):
# Line 627  class Defect(object): Line 627  class Defect(object):
627          :note: by default ``sqrt(self.bilinearform(x,x)`` is returned.          :note: by default ``sqrt(self.bilinearform(x,x)`` is returned.
628          """          """
629          s=self.bilinearform(x,x)          s=self.bilinearform(x,x)
630          if s<0: raise NegativeNorm,"negative norm."          if s<0: raise NegativeNorm("negative norm.")
631          return math.sqrt(s)          return math.sqrt(s)
632    
633      def eval(self,x):      def eval(self,x):
# Line 651  class Defect(object): Line 651  class Defect(object):
651          :param inc: relative increment length          :param inc: relative increment length
652          :type inc: positive ``float``          :type inc: positive ``float``
653          """          """
654          if inc<=0: raise ValueError,"positive increment required."          if inc<=0: raise ValueError("positive increment required.")
655          self.__inc=inc          self.__inc=inc
656    
657      def getDerivativeIncrementLength(self):      def getDerivativeIncrementLength(self):
# Line 722  def NewtonGMRES(defect, x, iter_max=100, Line 722  def NewtonGMRES(defect, x, iter_max=100,
722     :rtype: same type as the initial guess     :rtype: same type as the initial guess
723     """     """
724     lmaxit=iter_max     lmaxit=iter_max
725     if atol<0: raise ValueError,"atol needs to be non-negative."     if atol<0: raise ValueError("atol needs to be non-negative.")
726     if rtol<0: raise ValueError,"rtol needs to be non-negative."     if rtol<0: raise ValueError("rtol needs to be non-negative.")
727     if rtol+atol<=0: raise ValueError,"rtol or atol needs to be non-negative."     if rtol+atol<=0: raise ValueError("rtol or atol needs to be non-negative.")
728     if gamma<=0 or gamma>=1: raise ValueError,"tolerance safety factor for inner iteration (gamma =%s) needs to be positive and less than 1."%gamma     if gamma<=0 or gamma>=1: raise ValueError("tolerance safety factor for inner iteration (gamma =%s) needs to be positive and less than 1."%gamma)
729     if subtol_max<=0 or subtol_max>=1: raise ValueError,"upper bound for inner tolerance for inner iteration (subtol_max =%s) needs to be positive and less than 1."%subtol_max     if subtol_max<=0 or subtol_max>=1: raise ValueError("upper bound for inner tolerance for inner iteration (subtol_max =%s) needs to be positive and less than 1."%subtol_max)
730    
731     F=defect(x)     F=defect(x)
732     fnrm=defect.norm(F)     fnrm=defect.norm(F)
733     stop_tol=atol + rtol*fnrm     stop_tol=atol + rtol*fnrm
734     subtol=subtol_max     subtol=subtol_max
735     if verbose: print "NewtonGMRES: initial residual = %e."%fnrm     if verbose: print(("NewtonGMRES: initial residual = %e."%fnrm))
736     if verbose: print "             tolerance = %e."%subtol     if verbose: print(("             tolerance = %e."%subtol))
737     iter=1     iter=1
738     #     #
739     # main iteration loop     # main iteration loop
740     #     #
741     while not fnrm<=stop_tol:     while not fnrm<=stop_tol:
742              if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max              if iter  >= iter_max: raise MaxIterReached("maximum number of %s steps reached."%iter_max)
743              #              #
744          #   adjust subtol_          #   adjust subtol_
745          #          #
# Line 756  def NewtonGMRES(defect, x, iter_max=100, Line 756  def NewtonGMRES(defect, x, iter_max=100,
756              #     if  atol is reached sub_iter returns the numer of steps performed to get there              #     if  atol is reached sub_iter returns the numer of steps performed to get there
757              #              #
758              #              #
759              if verbose: print "             subiteration (GMRES) is called with relative tolerance %e."%subtol              if verbose: print(("             subiteration (GMRES) is called with relative tolerance %e."%subtol))
760              try:              try:
761                 xc, sub_iter=__FDGMRES(F, defect, x, subtol*fnrm, iter_max=iter_max-iter, iter_restart=sub_iter_max)                 xc, sub_iter=__FDGMRES(F, defect, x, subtol*fnrm, iter_max=iter_max-iter, iter_restart=sub_iter_max)
762              except MaxIterReached:              except MaxIterReached:
763                 raise MaxIterReached,"maximum number of %s steps reached."%iter_max                 raise MaxIterReached("maximum number of %s steps reached."%iter_max)
764              if sub_iter<0:              if sub_iter<0:
765                 iter+=sub_iter_max                 iter+=sub_iter_max
766              else:              else:
# Line 770  def NewtonGMRES(defect, x, iter_max=100, Line 770  def NewtonGMRES(defect, x, iter_max=100,
770              F=defect(x)              F=defect(x)
771          iter+=1          iter+=1
772              fnrmo, fnrm=fnrm, defect.norm(F)              fnrmo, fnrm=fnrm, defect.norm(F)
773              if verbose: print "             step %s: residual %e."%(iter,fnrm)              if verbose: print(("             step %s: residual %e."%(iter,fnrm)))
774     if verbose: print "NewtonGMRES: completed after %s steps."%iter     if verbose: print(("NewtonGMRES: completed after %s steps."%iter))
775     return x     return x
776    
777  def __givapp(c,s,vin):  def __givapp(c,s,vin):
# Line 807  def __FDGMRES(F0, defect, x0, atol, iter Line 807  def __FDGMRES(F0, defect, x0, atol, iter
807     iter=0     iter=0
808     while rho > atol and iter<iter_restart-1:     while rho > atol and iter<iter_restart-1:
809          if iter  >= iter_max:          if iter  >= iter_max:
810              raise MaxIterReached,"maximum number of %s steps reached."%iter_max              raise MaxIterReached("maximum number of %s steps reached."%iter_max)
811    
812          p=defect.derivative(F0,x0,v[iter], v_is_normalised=True)          p=defect.derivative(F0,x0,v[iter], v_is_normalised=True)
813          v.append(p)          v.append(p)
# Line 928  def GMRES(r, Aprod, x, bilinearform, ato Line 928  def GMRES(r, Aprod, x, bilinearform, ato
928     iter=0     iter=0
929     if rtol>0:     if rtol>0:
930        r_dot_r = bilinearform(r, r)        r_dot_r = bilinearform(r, r)
931        if r_dot_r<0: raise NegativeNorm,"negative norm."        if r_dot_r<0: raise NegativeNorm("negative norm.")
932        atol2=atol+rtol*math.sqrt(r_dot_r)        atol2=atol+rtol*math.sqrt(r_dot_r)
933        if verbose: print "GMRES: norm of right hand side = %e (absolute tolerance = %e)"%(math.sqrt(r_dot_r), atol2)        if verbose: print(("GMRES: norm of right hand side = %e (absolute tolerance = %e)"%(math.sqrt(r_dot_r), atol2)))
934     else:     else:
935        atol2=atol        atol2=atol
936        if verbose: print "GMRES: absolute tolerance = %e"%atol2        if verbose: print(("GMRES: absolute tolerance = %e"%atol2))
937     if atol2<=0:     if atol2<=0:
938        raise ValueError,"Non-positive tolarance."        raise ValueError("Non-positive tolarance.")
939    
940     while True:     while True:
941        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached"%iter_max        if iter  >= iter_max: raise MaxIterReached("maximum number of %s steps reached"%iter_max)
942        if restarted:        if restarted:
943           r2 = r-Aprod(x-x2)           r2 = r-Aprod(x-x2)
944        else:        else:
# Line 947  def GMRES(r, Aprod, x, bilinearform, ato Line 947  def GMRES(r, Aprod, x, bilinearform, ato
947        x,stopped=_GMRESm(r2, Aprod, x, bilinearform, atol2, iter_max=iter_max-iter, iter_restart=m, verbose=verbose,P_R=P_R)        x,stopped=_GMRESm(r2, Aprod, x, bilinearform, atol2, iter_max=iter_max-iter, iter_restart=m, verbose=verbose,P_R=P_R)
948        iter+=iter_restart        iter+=iter_restart
949        if stopped: break        if stopped: break
950        if verbose: print "GMRES: restart."        if verbose: print("GMRES: restart.")
951        restarted=True        restarted=True
952     if verbose: print "GMRES: tolerance has been reached."     if verbose: print("GMRES: tolerance has been reached.")
953     return x     return x
954    
955  def _GMRESm(r, Aprod, x, bilinearform, atol, iter_max=100, iter_restart=20, verbose=False, P_R=None):  def _GMRESm(r, Aprod, x, bilinearform, atol, iter_max=100, iter_restart=20, verbose=False, P_R=None):
# Line 962  def _GMRESm(r, Aprod, x, bilinearform, a Line 962  def _GMRESm(r, Aprod, x, bilinearform, a
962     v=[]     v=[]
963    
964     r_dot_r = bilinearform(r, r)     r_dot_r = bilinearform(r, r)
965     if r_dot_r<0: raise NegativeNorm,"negative norm."     if r_dot_r<0: raise NegativeNorm("negative norm.")
966     rho=math.sqrt(r_dot_r)     rho=math.sqrt(r_dot_r)
967    
968     v.append(r/rho)     v.append(r/rho)
969     g[0]=rho     g[0]=rho
970    
971     if verbose: print "GMRES: initial residual %e (absolute tolerance = %e)"%(rho,atol)     if verbose: print(("GMRES: initial residual %e (absolute tolerance = %e)"%(rho,atol)))
972     while not (rho<=atol or iter==iter_restart):     while not (rho<=atol or iter==iter_restart):
973    
974      if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max      if iter  >= iter_max: raise MaxIterReached("maximum number of %s steps reached."%iter_max)
975    
976          if P_R!=None:          if P_R!=None:
977              p=Aprod(P_R(v[iter]))              p=Aprod(P_R(v[iter]))
# Line 1018  def _GMRESm(r, Aprod, x, bilinearform, a Line 1018  def _GMRESm(r, Aprod, x, bilinearform, a
1018  # Update the residual norm  # Update the residual norm
1019    
1020          rho=abs(g[iter+1])          rho=abs(g[iter+1])
1021          if verbose: print "GMRES: iteration step %s: residual %e"%(iter,rho)          if verbose: print(("GMRES: iteration step %s: residual %e"%(iter,rho)))
1022      iter+=1      iter+=1
1023    
1024  # At this point either iter > iter_max or rho < tol.  # At this point either iter > iter_max or rho < tol.
1025  # It's time to compute x and leave.  # It's time to compute x and leave.
1026    
1027     if verbose: print "GMRES: iteration stopped after %s step."%iter     if verbose: print(("GMRES: iteration stopped after %s step."%iter))
1028     if iter > 0 :     if iter > 0 :
1029       y=numpy.zeros(iter,numpy.float64)       y=numpy.zeros(iter,numpy.float64)
1030       y[iter-1] = g[iter-1] / h[iter-1,iter-1]       y[iter-1] = g[iter-1] / h[iter-1,iter-1]
# Line 1108  def MINRES(r, Aprod, x, Msolve, bilinear Line 1108  def MINRES(r, Aprod, x, Msolve, bilinear
1108      y = Msolve(r)      y = Msolve(r)
1109      beta1 = bilinearform(y,r)      beta1 = bilinearform(y,r)
1110    
1111      if beta1< 0: raise NegativeNorm,"negative norm."      if beta1< 0: raise NegativeNorm("negative norm.")
1112    
1113      #  If r = 0 exactly, stop with x      #  If r = 0 exactly, stop with x
1114      if beta1==0: return x      if beta1==0: return x
# Line 1143  def MINRES(r, Aprod, x, Msolve, bilinear Line 1143  def MINRES(r, Aprod, x, Msolve, bilinear
1143      # --------------------------------------------------------------------      # --------------------------------------------------------------------
1144      while not rnorm<=atol+rtol*Anorm*ynorm:    #  checks ||r|| < (||A|| ||x||) * TOL      while not rnorm<=atol+rtol*Anorm*ynorm:    #  checks ||r|| < (||A|| ||x||) * TOL
1145    
1146      if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max      if iter  >= iter_max: raise MaxIterReached("maximum number of %s steps reached."%iter_max)
1147          iter    = iter  +  1          iter    = iter  +  1
1148    
1149          #-----------------------------------------------------------------          #-----------------------------------------------------------------
# Line 1173  def MINRES(r, Aprod, x, Msolve, bilinear Line 1173  def MINRES(r, Aprod, x, Msolve, bilinear
1173          y = Msolve(r2)          y = Msolve(r2)
1174          oldb   = beta                           # oldb = betak          oldb   = beta                           # oldb = betak
1175          beta   = bilinearform(y,r2)             # beta = betak+1^2          beta   = bilinearform(y,r2)             # beta = betak+1^2
1176          if beta < 0: raise NegativeNorm,"negative norm."          if beta < 0: raise NegativeNorm("negative norm.")
1177    
1178          beta   = math.sqrt( beta )          beta   = math.sqrt( beta )
1179          tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta          tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta
# Line 1280  def TFQMR(r, Aprod, x, bilinearform, ato Line 1280  def TFQMR(r, Aprod, x, bilinearform, ato
1280    theta = 0.0;    theta = 0.0;
1281    eta = 0.0;    eta = 0.0;
1282    rho=bilinearform(r,r)    rho=bilinearform(r,r)
1283    if rho < 0: raise NegativeNorm,"negative norm."    if rho < 0: raise NegativeNorm("negative norm.")
1284    tau = math.sqrt(rho)    tau = math.sqrt(rho)
1285    norm_r0=tau    norm_r0=tau
1286    while tau>atol+rtol*norm_r0:    while tau>atol+rtol*norm_r0:
1287      if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max      if iter  >= iter_max: raise MaxIterReached("maximum number of %s steps reached."%iter_max)
1288    
1289      sigma = bilinearform(r,v)      sigma = bilinearform(r,v)
1290      if sigma == 0.0: raise IterationBreakDown,'TFQMR breakdown, sigma=0'      if sigma == 0.0: raise IterationBreakDown('TFQMR breakdown, sigma=0')
1291    
1292      alpha = rho / sigma      alpha = rho / sigma
1293    
# Line 1315  def TFQMR(r, Aprod, x, bilinearform, ato Line 1315  def TFQMR(r, Aprod, x, bilinearform, ato
1315  #  #
1316  #  Try to terminate the iteration at each pass through the loop  #  Try to terminate the iteration at each pass through the loop
1317  #  #
1318      if rho == 0.0: raise IterationBreakDown,'TFQMR breakdown, rho=0'      if rho == 0.0: raise IterationBreakDown('TFQMR breakdown, rho=0')
1319    
1320      rhon = bilinearform(r,w)      rhon = bilinearform(r,w)
1321      beta = rhon / rho;      beta = rhon / rho;
# Line 1387  class ArithmeticTuple(object): Line 1387  class ArithmeticTuple(object):
1387         try:         try:
1388             l=len(other)             l=len(other)
1389             if l!=len(self):             if l!=len(self):
1390                 raise ValueError,"length of arguments don't match."                 raise ValueError("length of arguments don't match.")
1391             for i in range(l): out.append(self[i]*other[i])             for i in range(l): out.append(self[i]*other[i])
1392         except TypeError:         except TypeError:
1393             for i in range(len(self)): out.append(self[i]*other)             for i in range(len(self)): out.append(self[i]*other)
# Line 1406  class ArithmeticTuple(object): Line 1406  class ArithmeticTuple(object):
1406         try:         try:
1407             l=len(other)             l=len(other)
1408             if l!=len(self):             if l!=len(self):
1409                 raise ValueError,"length of arguments don't match."                 raise ValueError("length of arguments don't match.")
1410             for i in range(l): out.append(other[i]*self[i])             for i in range(l): out.append(other[i]*self[i])
1411         except TypeError:         except TypeError:
1412             for i in range(len(self)): out.append(other*self[i])             for i in range(len(self)): out.append(other*self[i])
# Line 1436  class ArithmeticTuple(object): Line 1436  class ArithmeticTuple(object):
1436         try:         try:
1437             l=len(other)             l=len(other)
1438             if l!=len(self):             if l!=len(self):
1439                 raise ValueError,"length of arguments don't match."                 raise ValueError("length of arguments don't match.")
1440             for i in range(l): out.append(other[i]/self[i])             for i in range(l): out.append(other[i]/self[i])
1441         except TypeError:         except TypeError:
1442             for i in range(len(self)): out.append(other/self[i])             for i in range(len(self)): out.append(other/self[i])
# Line 1450  class ArithmeticTuple(object): Line 1450  class ArithmeticTuple(object):
1450         :type other: ``ArithmeticTuple``         :type other: ``ArithmeticTuple``
1451         """         """
1452         if len(self) != len(other):         if len(self) != len(other):
1453             raise ValueError,"tuple lengths must match."             raise ValueError("tuple lengths must match.")
1454         for i in range(len(self)):         for i in range(len(self)):
1455             self.__items[i]+=other[i]             self.__items[i]+=other[i]
1456         return self         return self
# Line 1466  class ArithmeticTuple(object): Line 1466  class ArithmeticTuple(object):
1466         try:         try:
1467             l=len(other)             l=len(other)
1468             if l!=len(self):             if l!=len(self):
1469                 raise ValueError,"length of arguments don't match."                 raise ValueError("length of arguments don't match.")
1470             for i in range(l): out.append(self[i]+other[i])             for i in range(l): out.append(self[i]+other[i])
1471         except TypeError:         except TypeError:
1472             for i in range(len(self)): out.append(self[i]+other)             for i in range(len(self)): out.append(self[i]+other)
# Line 1483  class ArithmeticTuple(object): Line 1483  class ArithmeticTuple(object):
1483         try:         try:
1484             l=len(other)             l=len(other)
1485             if l!=len(self):             if l!=len(self):
1486                 raise ValueError,"length of arguments don't match."                 raise ValueError("length of arguments don't match.")
1487             for i in range(l): out.append(self[i]-other[i])             for i in range(l): out.append(self[i]-other[i])
1488         except TypeError:         except TypeError:
1489             for i in range(len(self)): out.append(self[i]-other)             for i in range(len(self)): out.append(self[i]-other)
# Line 1497  class ArithmeticTuple(object): Line 1497  class ArithmeticTuple(object):
1497         :type other: ``ArithmeticTuple``         :type other: ``ArithmeticTuple``
1498         """         """
1499         if len(self) != len(other):         if len(self) != len(other):
1500             raise ValueError,"tuple length must match."             raise ValueError("tuple length must match.")
1501         for i in range(len(self)):         for i in range(len(self)):
1502             self.__items[i]-=other[i]             self.__items[i]-=other[i]
1503         return self         return self
# Line 1565  class HomogeneousSaddlePointProblem(obje Line 1565  class HomogeneousSaddlePointProblem(obje
1565           """           """
1566           if not K_p == None:           if not K_p == None:
1567              if K_p<1:              if K_p<1:
1568                 raise ValueError,"K_p need to be greater or equal to 1."                 raise ValueError("K_p need to be greater or equal to 1.")
1569           else:           else:
1570              K_p=self.__K_p              K_p=self.__K_p
1571    
1572           if not K_v == None:           if not K_v == None:
1573              if K_v<1:              if K_v<1:
1574                 raise ValueError,"K_v need to be greater or equal to 1."                 raise ValueError("K_v need to be greater or equal to 1.")
1575           else:           else:
1576              K_v=self.__K_v              K_v=self.__K_v
1577    
1578           if not rtol_max == None:           if not rtol_max == None:
1579              if rtol_max<=0 or rtol_max>=1:              if rtol_max<=0 or rtol_max>=1:
1580                 raise ValueError,"rtol_max needs to be positive and less than 1."                 raise ValueError("rtol_max needs to be positive and less than 1.")
1581           else:           else:
1582              rtol_max=self.__rtol_max              rtol_max=self.__rtol_max
1583    
1584           if not rtol_min == None:           if not rtol_min == None:
1585              if rtol_min<=0 or rtol_min>=1:              if rtol_min<=0 or rtol_min>=1:
1586                 raise ValueError,"rtol_min needs to be positive and less than 1."                 raise ValueError("rtol_min needs to be positive and less than 1.")
1587           else:           else:
1588              rtol_min=self.__rtol_min              rtol_min=self.__rtol_min
1589    
1590           if not chi_max == None:           if not chi_max == None:
1591              if chi_max<=0 or chi_max>=1:              if chi_max<=0 or chi_max>=1:
1592                 raise ValueError,"chi_max needs to be positive and less than 1."                 raise ValueError("chi_max needs to be positive and less than 1.")
1593           else:           else:
1594              chi_max = self.__chi_max              chi_max = self.__chi_max
1595    
1596           if not reduction_factor == None:           if not reduction_factor == None:
1597              if reduction_factor<=0 or reduction_factor>1:              if reduction_factor<=0 or reduction_factor>1:
1598                 raise ValueError,"reduction_factor need to be between zero and one."                 raise ValueError("reduction_factor need to be between zero and one.")
1599           else:           else:
1600              reduction_factor=self.__reduction_factor              reduction_factor=self.__reduction_factor
1601    
1602           if not theta == None:           if not theta == None:
1603              if theta<=0 or theta>1:              if theta<=0 or theta>1:
1604                 raise ValueError,"theta need to be between zero and one."                 raise ValueError("theta need to be between zero and one.")
1605           else:           else:
1606              theta=self.__theta              theta=self.__theta
1607    
1608           if rtol_min>=rtol_max:           if rtol_min>=rtol_max:
1609               raise ValueError,"rtol_max = %e needs to be greater than rtol_min = %e"%(rtol_max,rtol_min)               raise ValueError("rtol_max = %e needs to be greater than rtol_min = %e"%(rtol_max,rtol_min))
1610           self.__chi_max = chi_max           self.__chi_max = chi_max
1611           self.__rtol_max = rtol_max           self.__rtol_max = rtol_max
1612           self.__K_p = K_p           self.__K_p = K_p
# Line 1626  class HomogeneousSaddlePointProblem(obje Line 1626  class HomogeneousSaddlePointProblem(obje
1626           :rtype: ``float``           :rtype: ``float``
1627           :note: used if PCG is applied.           :note: used if PCG is applied.
1628           """           """
1629           raise NotImplementedError,"no inner product for p and Bv implemented."           raise NotImplementedError("no inner product for p and Bv implemented.")
1630    
1631        def inner_p(self,p0,p1):        def inner_p(self,p0,p1):
1632           """           """
# Line 1637  class HomogeneousSaddlePointProblem(obje Line 1637  class HomogeneousSaddlePointProblem(obje
1637           :return: inner product of p0 and p1           :return: inner product of p0 and p1
1638           :rtype: ``float``           :rtype: ``float``
1639           """           """
1640           raise NotImplementedError,"no inner product for p implemented."           raise NotImplementedError("no inner product for p implemented.")
1641        
1642        def norm_v(self,v):        def norm_v(self,v):
1643           """           """
# Line 1647  class HomogeneousSaddlePointProblem(obje Line 1647  class HomogeneousSaddlePointProblem(obje
1647           :return: norm of v           :return: norm of v
1648           :rtype: non-negative ``float``           :rtype: non-negative ``float``
1649           """           """
1650           raise NotImplementedError,"no norm of v implemented."           raise NotImplementedError("no norm of v implemented.")
1651        def getDV(self, p, v, tol):        def getDV(self, p, v, tol):
1652           """           """
1653           return a correction to the value for a given v and a given p with accuracy `tol` (overwrite)           return a correction to the value for a given v and a given p with accuracy `tol` (overwrite)
# Line 1657  class HomogeneousSaddlePointProblem(obje Line 1657  class HomogeneousSaddlePointProblem(obje
1657           :return: dv given as *dv= A^{-1} (f-A v-B^*p)*           :return: dv given as *dv= A^{-1} (f-A v-B^*p)*
1658           :note: Only *A* may depend on *v* and *p*           :note: Only *A* may depend on *v* and *p*
1659           """           """
1660           raise NotImplementedError,"no dv calculation implemented."           raise NotImplementedError("no dv calculation implemented.")
1661    
1662                    
1663        def Bv(self,v, tol):        def Bv(self,v, tol):
# Line 1667  class HomogeneousSaddlePointProblem(obje Line 1667  class HomogeneousSaddlePointProblem(obje
1667          :rtype: equal to the type of p          :rtype: equal to the type of p
1668          :note: boundary conditions on p should be zero!          :note: boundary conditions on p should be zero!
1669          """          """
1670          raise NotImplementedError, "no operator B implemented."          raise NotImplementedError("no operator B implemented.")
1671    
1672        def norm_Bv(self,Bv):        def norm_Bv(self,Bv):
1673          """          """
# Line 1676  class HomogeneousSaddlePointProblem(obje Line 1676  class HomogeneousSaddlePointProblem(obje
1676          :rtype: equal to the type of p          :rtype: equal to the type of p
1677          :note: boundary conditions on p should be zero!          :note: boundary conditions on p should be zero!
1678          """          """
1679          raise NotImplementedError, "no norm of Bv implemented."          raise NotImplementedError("no norm of Bv implemented.")
1680    
1681        def solve_AinvBt(self,dp, tol):        def solve_AinvBt(self,dp, tol):
1682           """           """
# Line 1686  class HomogeneousSaddlePointProblem(obje Line 1686  class HomogeneousSaddlePointProblem(obje
1686           :return: the solution of *A dv=B^*dp*           :return: the solution of *A dv=B^*dp*
1687           :note: boundary conditions on dv should be zero! *A* is the operator used in ``getDV`` and must not be altered.           :note: boundary conditions on dv should be zero! *A* is the operator used in ``getDV`` and must not be altered.
1688           """           """
1689           raise NotImplementedError,"no operator A implemented."           raise NotImplementedError("no operator A implemented.")
1690    
1691        def solve_prec(self,Bv, tol):        def solve_prec(self,Bv, tol):
1692           """           """
# Line 1695  class HomogeneousSaddlePointProblem(obje Line 1695  class HomogeneousSaddlePointProblem(obje
1695           :rtype: equal to the type of p           :rtype: equal to the type of p
1696           :note: boundary conditions on p should be zero!           :note: boundary conditions on p should be zero!
1697           """           """
1698           raise NotImplementedError,"no preconditioner for Schur complement implemented."           raise NotImplementedError("no preconditioner for Schur complement implemented.")
1699        #=============================================================        #=============================================================
1700        def __Aprod_PCG(self,dp):        def __Aprod_PCG(self,dp):
1701            dv=self.solve_AinvBt(dp, self.__subtol)            dv=self.solve_AinvBt(dp, self.__subtol)
# Line 1722  class HomogeneousSaddlePointProblem(obje Line 1722  class HomogeneousSaddlePointProblem(obje
1722            :rtype: ``float``            :rtype: ``float``
1723            """            """
1724            f=self.inner_p(p,p)            f=self.inner_p(p,p)
1725            if f<0: raise ValueError,"negative pressure norm."            if f<0: raise ValueError("negative pressure norm.")
1726            return math.sqrt(f)            return math.sqrt(f)
1727                
1728        def solve(self,v,p,max_iter=20, verbose=False, usePCG=True, iter_restart=20, max_correction_steps=10):        def solve(self,v,p,max_iter=20, verbose=False, usePCG=True, iter_restart=20, max_correction_steps=10):
# Line 1764  class HomogeneousSaddlePointProblem(obje Line 1764  class HomogeneousSaddlePointProblem(obje
1764           chi=None           chi=None
1765           eps=None           eps=None
1766    
1767           if self.verbose: print "HomogeneousSaddlePointProblem: start iteration: rtol= %e, atol=%e"%(rtol, atol)           if self.verbose: print(("HomogeneousSaddlePointProblem: start iteration: rtol= %e, atol=%e"%(rtol, atol)))
1768           while not converged:           while not converged:
1769    
1770               # get tolerance for velecity increment:               # get tolerance for velecity increment:
# Line 1773  class HomogeneousSaddlePointProblem(obje Line 1773  class HomogeneousSaddlePointProblem(obje
1773               else:               else:
1774                  rtol_v=min(chi/K_v,self.__rtol_max)                  rtol_v=min(chi/K_v,self.__rtol_max)
1775               rtol_v=max(rtol_v, self.__rtol_min)               rtol_v=max(rtol_v, self.__rtol_min)
1776               if self.verbose: print "HomogeneousSaddlePointProblem: step %s: rtol_v= %e"%(correction_step,rtol_v)               if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: rtol_v= %e"%(correction_step,rtol_v)))
1777               # get velocity increment:               # get velocity increment:
1778               dv1=self.getDV(p,v,rtol_v)               dv1=self.getDV(p,v,rtol_v)
1779               v1=v+dv1               v1=v+dv1
1780               Bv1=self.Bv(v1, rtol_v)               Bv1=self.Bv(v1, rtol_v)
1781               norm_Bv1=self.norm_Bv(Bv1)               norm_Bv1=self.norm_Bv(Bv1)
1782               norm_dv1=self.norm_v(dv1)               norm_dv1=self.norm_v(dv1)
1783               if self.verbose: print "HomogeneousSaddlePointProblem: step %s: norm_Bv1 = %e, norm_dv1 = %e"%(correction_step, norm_Bv1, norm_dv1)               if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: norm_Bv1 = %e, norm_dv1 = %e"%(correction_step, norm_Bv1, norm_dv1)))
1784               if norm_dv1*self.__theta < norm_Bv1:               if norm_dv1*self.__theta < norm_Bv1:
1785                  # get tolerance for pressure increment:                  # get tolerance for pressure increment:
1786                  large_Bv1=True                  large_Bv1=True
# Line 1789  class HomogeneousSaddlePointProblem(obje Line 1789  class HomogeneousSaddlePointProblem(obje
1789                  else:                  else:
1790                     rtol_p=min(chi**2*eps/K_p/norm_Bv1, self.__rtol_max)                     rtol_p=min(chi**2*eps/K_p/norm_Bv1, self.__rtol_max)
1791                  self.__subtol=max(rtol_p**2, self.__rtol_min)                  self.__subtol=max(rtol_p**2, self.__rtol_min)
1792                  if self.verbose: print "HomogeneousSaddlePointProblem: step %s: rtol_p= %e"%(correction_step,rtol_p)                  if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: rtol_p= %e"%(correction_step,rtol_p)))
1793                  # now we solve for the pressure increment dp from B*A^{-1}B^* dp = Bv1                  # now we solve for the pressure increment dp from B*A^{-1}B^* dp = Bv1
1794                  if usePCG:                  if usePCG:
1795                      dp,r,a_norm=PCG(ArithmeticTuple(v1,Bv1),self.__Aprod_PCG,0*p,self.__Msolve_PCG,self.__inner_PCG,atol=0, rtol=rtol_p,iter_max=max_iter, verbose=self.verbose)                      dp,r,a_norm=PCG(ArithmeticTuple(v1,Bv1),self.__Aprod_PCG,0*p,self.__Msolve_PCG,self.__inner_PCG,atol=0, rtol=rtol_p,iter_max=max_iter, verbose=self.verbose)
# Line 1809  class HomogeneousSaddlePointProblem(obje Line 1809  class HomogeneousSaddlePointProblem(obje
1809               # update business:               # update business:
1810               norm_dv2=self.norm_v(v2-v)               norm_dv2=self.norm_v(v2-v)
1811               norm_v2=self.norm_v(v2)               norm_v2=self.norm_v(v2)
1812               if self.verbose: print "HomogeneousSaddlePointProblem: step %s: v2 = %e, norm_dv2 = %e"%(correction_step, norm_v2, self.norm_v(v2-v))               if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: v2 = %e, norm_dv2 = %e"%(correction_step, norm_v2, self.norm_v(v2-v))))
1813               eps, eps_old = max(norm_Bv1, norm_dv2), eps               eps, eps_old = max(norm_Bv1, norm_dv2), eps
1814               if eps_old == None:               if eps_old == None:
1815                    chi, chi_old = None, chi                    chi, chi_old = None, chi
# Line 1817  class HomogeneousSaddlePointProblem(obje Line 1817  class HomogeneousSaddlePointProblem(obje
1817                    chi, chi_old = min(eps/ eps_old, self.__chi_max), chi                    chi, chi_old = min(eps/ eps_old, self.__chi_max), chi
1818               if eps != None:               if eps != None:
1819                   if chi !=None:                   if chi !=None:
1820                      if self.verbose: print "HomogeneousSaddlePointProblem: step %s: convergence rate = %e, correction = %e"%(correction_step,chi, eps)                      if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: convergence rate = %e, correction = %e"%(correction_step,chi, eps)))
1821                   else:                   else:
1822                      if self.verbose: print "HomogeneousSaddlePointProblem: step %s: correction = %e"%(correction_step, eps)                      if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: correction = %e"%(correction_step, eps)))
1823               if eps <= rtol*norm_v2+atol :               if eps <= rtol*norm_v2+atol :
1824                   converged = True                   converged = True
1825               else:               else:
1826                   if correction_step>=max_correction_steps:                   if correction_step>=max_correction_steps:
1827                        raise CorrectionFailed,"Given up after %d correction steps."%correction_step                        raise CorrectionFailed("Given up after %d correction steps."%correction_step)
1828                   if chi_old!=None:                   if chi_old!=None:
1829                      K_p=max(1,self.__reduction_factor*K_p,(chi-chi_old)/chi_old**2*K_p)                      K_p=max(1,self.__reduction_factor*K_p,(chi-chi_old)/chi_old**2*K_p)
1830                      K_v=max(1,self.__reduction_factor*K_v,(chi-chi_old)/chi_old**2*K_p)                      K_v=max(1,self.__reduction_factor*K_v,(chi-chi_old)/chi_old**2*K_p)
1831                      if self.verbose: print "HomogeneousSaddlePointProblem: step %s: new adjustment factor K = %e"%(correction_step,K_p)                      if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: new adjustment factor K = %e"%(correction_step,K_p)))
1832               correction_step+=1               correction_step+=1
1833               v,p =v2, p2               v,p =v2, p2
1834           if self.verbose: print "HomogeneousSaddlePointProblem: tolerance reached after %s steps."%correction_step           if self.verbose: print(("HomogeneousSaddlePointProblem: tolerance reached after %s steps."%correction_step))
1835       return v,p       return v,p
1836        #========================================================================        #========================================================================
1837        def setTolerance(self,tolerance=1.e-4):        def setTolerance(self,tolerance=1.e-4):
# Line 1842  class HomogeneousSaddlePointProblem(obje Line 1842  class HomogeneousSaddlePointProblem(obje
1842           :type tolerance: non-negative ``float``           :type tolerance: non-negative ``float``
1843           """           """
1844           if tolerance<0:           if tolerance<0:
1845               raise ValueError,"tolerance must be positive."               raise ValueError("tolerance must be positive.")
1846           self.__rtol=tolerance           self.__rtol=tolerance
1847    
1848        def getTolerance(self):        def getTolerance(self):
# Line 1862  class HomogeneousSaddlePointProblem(obje Line 1862  class HomogeneousSaddlePointProblem(obje
1862           :type tolerance: non-negative ``float``           :type tolerance: non-negative ``float``
1863           """           """
1864           if tolerance<0:           if tolerance<0:
1865               raise ValueError,"tolerance must be non-negative."               raise ValueError("tolerance must be non-negative.")
1866           self.__atol=tolerance           self.__atol=tolerance
1867    
1868        def getAbsoluteTolerance(self):        def getAbsoluteTolerance(self):

Legend:
Removed from v.3676  
changed lines
  Added in v.3771

  ViewVC Help
Powered by ViewVC 1.1.26