/[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 2793 by gross, Tue Dec 1 06:10:10 2009 UTC revision 3911 by jfenwick, Thu Jun 14 01:01:03 2012 UTC
# Line 1  Line 1 
1    
2  ########################################################  ########################################################
3  #  #
4  # Copyright (c) 2003-2009 by University of Queensland  # Copyright (c) 2003-2012 by University of Queensland
5  # Earth Systems Science Computational Center (ESSCC)  # Earth Systems Science Computational Center (ESSCC)
6  # http://www.uq.edu.au/esscc  # http://www.uq.edu.au/esscc
7  #  #
# Line 11  Line 11 
11  #  #
12  ########################################################  ########################################################
13    
14  __copyright__="""Copyright (c) 2003-2009 by University of Queensland  __copyright__="""Copyright (c) 2003-2012 by University of Queensland
15  Earth Systems Science Computational Center (ESSCC)  Earth Systems Science Computational Center (ESSCC)
16  http://www.uq.edu.au/esscc  http://www.uq.edu.au/esscc
17  Primary Business: Queensland, Australia"""  Primary Business: Queensland, Australia"""
# Line 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 153  class Projector: Line 153  class Projector:
153      """      """
154      return self.__pde.getSolverOptions()      return self.__pde.getSolverOptions()
155    
156      def getValue(self, input_data):
157        """
158        Projects ``input_data`` onto a continuous function.
159    
160        :param input_data: the data to be projected
161        """
162        return self(input_data)
163    
164    def __call__(self, input_data):    def __call__(self, input_data):
165      """      """
166      Projects ``input_data`` onto a continuous function.      Projects ``input_data`` onto a continuous function.
# Line 291  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 332  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
347             except TypeError:             except TypeError:
348               iterative=False               iterative=False
349           xxx=self.__function_space.getX()
350         if iterative:         if iterative:
351             self.__id=[]             self.__id=[]
352             for p in x:             for p in x:
353                self.__id.append(util.length(self.__function_space.getX()-p[:self.__function_space.getDim()]).minGlobalDataPoint())                self.__id.append(util.length(xxx-p[:self.__function_space.getDim()]).minGlobalDataPoint())
354         else:         else:
355             self.__id=util.length(self.__function_space.getX()-x[:self.__function_space.getDim()]).minGlobalDataPoint()             self.__id=util.length(xxx-x[:self.__function_space.getDim()]).minGlobalDataPoint()
356    
357       def __str__(self):       def __str__(self):
358         """         """
# Line 421  class Locator: Line 430  class Locator:
430                  return out                  return out
431          else:          else:
432             return data             return data
433              
434    #     def setValue(self, data, v):
435    #       """
436    #       Sets the value of the ``data`` at the Locator.
437    #       """
438    #       data.expand()   # Need to ensure that this is done globally
439    #       if isinstance(data, escript.Data):
440    #    id=self.getId()
441    #    if isinstance(id, list):
442    #      for i in id:
443    #         data._setTupleForGlobalDataPoint(i[1], i[0], v)
444    #    else:
445    #      data._setTupleForGlobalDataPoint(id[1], id[0], v)
446    #       else:
447    #    raise TypeError, "setValue: Invalid argument type."
448    
449    
450  def getInfLocator(arg):  def getInfLocator(arg):
# Line 428  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 440  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 539  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)
582         x += alpha * d         x += alpha * d
583         if isinstance(q,ArithmeticTuple):         if isinstance(q,ArithmeticTuple):
584         r += q * (-alpha)      # Doing it the other way calls the float64.__mul__ not AT.__rmul__            r += q * (-alpha)      # Doing it the other way calls the float64.__mul__ not AT.__rmul__
585         else:         else:
586             r += (-alpha) * q             r += (-alpha) * q
587         rhat=Msolve(r)         rhat=Msolve(r)
# Line 567  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 603  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 627  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 698  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          #          #
746              if iter > 1:              if iter > 1:
747             rat=fnrm/fnrmo                 rat=fnrm/fnrmo
748                 subtol_old=subtol                 subtol_old=subtol
749             subtol=gamma*rat**2                 subtol=gamma*rat**2
750             if gamma*subtol_old**2 > .1: subtol=max(subtol,gamma*subtol_old**2)                 if gamma*subtol_old**2 > .1: subtol=max(subtol,gamma*subtol_old**2)
751             subtol=max(min(subtol,subtol_max), .5*stop_tol/fnrm)                 subtol=max(min(subtol,subtol_max), .5*stop_tol/fnrm)
752          #          #
753          # calculate newton increment xc          # calculate newton increment xc
754              #     if iter_max in __FDGMRES is reached MaxIterReached is thrown              #     if iter_max in __FDGMRES is reached MaxIterReached is thrown
# Line 732  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:
767                 iter+=sub_iter                 iter+=sub_iter
768              # ====              # ====
769          x+=xc              x+=xc
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 763  def __givapp(c,s,vin): Line 787  def __givapp(c,s,vin):
787      else:      else:
788          for i in range(len(c)):          for i in range(len(c)):
789              w1=c[i]*vrot[i]-s[i]*vrot[i+1]              w1=c[i]*vrot[i]-s[i]*vrot[i+1]
790          w2=s[i]*vrot[i]+c[i]*vrot[i+1]              w2=s[i]*vrot[i]+c[i]*vrot[i+1]
791              vrot[i]=w1              vrot[i]=w1
792              vrot[i+1]=w2              vrot[i+1]=w2
793      return vrot      return vrot
# Line 783  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 845  def __FDGMRES(F0, defect, x0, atol, iter Line 869  def __FDGMRES(F0, defect, x0, atol, iter
869            i=i-1            i=i-1
870       xhat=v[iter-1]*y[iter-1]       xhat=v[iter-1]*y[iter-1]
871       for i in range(iter-1):       for i in range(iter-1):
872      xhat += v[i]*y[i]         xhat += v[i]*y[i]
873     else :     else :
874        xhat=v[0] * 0        xhat=v[0] * 0
875    
# Line 904  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 923  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 938  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]))
978          else:          else:
979          p=Aprod(v[iter])              p=Aprod(v[iter])
980      v.append(p)          v.append(p)
981    
982      v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))          v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
983    
984  # Modified Gram-Schmidt  # Modified Gram-Schmidt
985      for j in range(iter+1):          for j in range(iter+1):
986        h[j,iter]=bilinearform(v[j],v[iter+1])            h[j,iter]=bilinearform(v[j],v[iter+1])
987        v[iter+1]-=h[j,iter]*v[j]            v[iter+1]-=h[j,iter]*v[j]
988    
989      h[iter+1,iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))          h[iter+1,iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
990      v_norm2=h[iter+1,iter]          v_norm2=h[iter+1,iter]
991    
992  # Reorthogonalize if needed  # Reorthogonalize if needed
993      if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)          if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
994       for j in range(iter+1):           for j in range(iter+1):
995          hr=bilinearform(v[j],v[iter+1])              hr=bilinearform(v[j],v[iter+1])
996              h[j,iter]=h[j,iter]+hr              h[j,iter]=h[j,iter]+hr
997              v[iter+1] -= hr*v[j]              v[iter+1] -= hr*v[j]
998    
999       v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))           v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
1000       h[iter+1,iter]=v_norm2           h[iter+1,iter]=v_norm2
1001    
1002  #   watch out for happy breakdown  #   watch out for happy breakdown
1003          if not v_norm2 == 0:          if not v_norm2 == 0:
1004           v[iter+1]=v[iter+1]/h[iter+1,iter]           v[iter+1]=v[iter+1]/h[iter+1,iter]
1005    
1006  #   Form and store the information for the new Givens rotation  #   Form and store the information for the new Givens rotation
1007      if iter > 0: h[:iter+1,iter]=__givapp(c[:iter],s[:iter],h[:iter+1,iter])          if iter > 0: h[:iter+1,iter]=__givapp(c[:iter],s[:iter],h[:iter+1,iter])
1008      mu=math.sqrt(h[iter,iter]*h[iter,iter]+h[iter+1,iter]*h[iter+1,iter])          mu=math.sqrt(h[iter,iter]*h[iter,iter]+h[iter+1,iter]*h[iter+1,iter])
1009    
1010      if mu!=0 :          if mu!=0 :
1011          c[iter]=h[iter,iter]/mu                  c[iter]=h[iter,iter]/mu
1012          s[iter]=-h[iter+1,iter]/mu                  s[iter]=-h[iter+1,iter]/mu
1013          h[iter,iter]=c[iter]*h[iter,iter]-s[iter]*h[iter+1,iter]                  h[iter,iter]=c[iter]*h[iter,iter]-s[iter]*h[iter+1,iter]
1014          h[iter+1,iter]=0.0                  h[iter+1,iter]=0.0
1015                  gg=__givapp(c[iter],s[iter],[g[iter],g[iter+1]])                  gg=__givapp(c[iter],s[iter],[g[iter],g[iter+1]])
1016                  g[iter]=gg[0]                  g[iter]=gg[0]
1017                  g[iter+1]=gg[1]                  g[iter+1]=gg[1]
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 1011  def _GMRESm(r, Aprod, x, bilinearform, a Line 1035  def _GMRESm(r, Aprod, x, bilinearform, a
1035            i=i-1            i=i-1
1036       xhat=v[iter-1]*y[iter-1]       xhat=v[iter-1]*y[iter-1]
1037       for i in range(iter-1):       for i in range(iter-1):
1038      xhat += v[i]*y[i]         xhat += v[i]*y[i]
1039     else:     else:
1040       xhat=v[0] * 0       xhat=v[0] * 0
1041     if P_R!=None:     if P_R!=None:
# Line 1084  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 1119  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 1149  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 1256  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 1291  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 1363  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 1382  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 1412  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 1426  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 1442  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 1459  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 1473  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 1501  class HomogeneousSaddlePointProblem(obje Line 1525  class HomogeneousSaddlePointProblem(obje
1525        *A* may depend weakly on *v* and *p*.        *A* may depend weakly on *v* and *p*.
1526        """        """
1527        def __init__(self, **kwargs):        def __init__(self, **kwargs):
1528      """          """
1529      initializes the saddle point problem          initializes the saddle point problem
1530      """          """
1531          self.resetControlParameters()          self.resetControlParameters()
1532          self.setTolerance()          self.setTolerance()
1533          self.setAbsoluteTolerance()          self.setAbsoluteTolerance()
1534        def resetControlParameters(self,gamma=0.85, gamma_min=1.e-2,chi_max=0.1, omega_div=0.2, omega_conv=1.1, rtol_min=1.e-7, rtol_max=0.9, chi=1., C_p=1., C_v=1., safety_factor=0.3):        def resetControlParameters(self, K_p=1., K_v=1., rtol_max=0.01, rtol_min = 1.e-7, chi_max=0.5, reduction_factor=0.3, theta = 0.1):
1535           """           """
1536           sets a control parameter           sets a control parameter
1537    
1538           :param gamma: ``1/(1-gamma)`` controls the perturbation of the converegence rate due to termination errors in the subproblems.           :param K_p: initial value for constant to adjust pressure tolerance
1539           :type gamma: ``float``           :type K_p: ``float``
1540           :param gamma_min: minimum value for ``gamma``.           :param K_v: initial value for constant to adjust velocity tolerance
1541           :type gamma_min: ``float``           :type K_v: ``float``
          :param chi_max: maximum tolerable converegence rate.  
          :type chi_max: ``float``  
          :param omega_div: reduction fact for ``gamma`` if no convergence is detected.  
          :type omega_div: ``float``  
          :param omega_conv: raise fact for ``gamma`` if convergence is detected.  
          :type omega_conv: ``float``  
          :param rtol_min: minimum relative tolerance used to calculate presssure and velocity increment.  
          :type rtol_min: ``float``  
1542           :param rtol_max: maximuim relative tolerance used to calculate presssure and velocity increment.           :param rtol_max: maximuim relative tolerance used to calculate presssure and velocity increment.
1543           :type rtol_max: ``float``           :type rtol_max: ``float``
1544           :param chi: initial convergence measure.           :param chi_max: maximum tolerable converegence rate.
1545           :type chi: ``float``           :type chi_max: ``float``
1546           :param C_p: initial value for constant to adjust pressure tolerance           :param reduction_factor: reduction factor for adjustment factors.
1547           :type C_p: ``float``           :type reduction_factor: ``float``
          :param C_v: initial value for constant to adjust velocity tolerance  
          :type C_v: ``float``  
          :param safety_factor: safety factor for addjustment of pressure and velocity tolerance from stopping criteria  
          :type safety_factor: ``float``  
1548           """           """
1549           self.setControlParameter(gamma, gamma_min ,chi_max , omega_div , omega_conv, rtol_min , rtol_max, chi,C_p, C_v,safety_factor)           self.setControlParameter(K_p, K_v, rtol_max, rtol_min, chi_max, reduction_factor, theta)
1550    
1551        def setControlParameter(self,gamma=None, gamma_min=None ,chi_max=None, omega_div=None, omega_conv=None, rtol_min=None, rtol_max=None, chi=None, C_p=None, C_v=None, safety_factor=None):        def setControlParameter(self,K_p=None, K_v=None, rtol_max=None, rtol_min=None, chi_max=None, reduction_factor=None, theta=None):
1552           """           """
1553           sets a control parameter           sets a control parameter
1554    
1555           :param gamma: ``1/(1-gamma)`` controls the perturbation of the converegence rate due to termination errors in the subproblems.  
1556           :type gamma: ``float``           :param K_p: initial value for constant to adjust pressure tolerance
1557           :param gamma_min: minimum value for ``gamma``.           :type K_p: ``float``
1558           :type gamma_min: ``float``           :param K_v: initial value for constant to adjust velocity tolerance
1559           :param chi_max: maximum tolerable converegence rate.           :type K_v: ``float``
          :type chi_max: ``float``  
          :param omega_div: reduction fact for ``gamma`` if no convergence is detected.  
          :type omega_div: ``float``  
          :param omega_conv: raise fact for ``gamma`` if convergence is detected.  
          :type omega_conv: ``float``  
          :param rtol_min: minimum relative tolerance used to calculate presssure and velocity increment.  
          :type rtol_min: ``float``  
1560           :param rtol_max: maximuim relative tolerance used to calculate presssure and velocity increment.           :param rtol_max: maximuim relative tolerance used to calculate presssure and velocity increment.
1561           :type rtol_max: ``float``           :type rtol_max: ``float``
1562           :param chi: initial convergence measure.           :param chi_max: maximum tolerable converegence rate.
1563           :type chi: ``float``           :type chi_max: ``float``
1564           :param C_p: initial value for constant to adjust pressure tolerance           :type reduction_factor: ``float``
1565           :type C_p: ``float``           """
1566           :param C_v: initial value for constant to adjust velocity tolerance           if not K_p == None:
1567           :type C_v: ``float``              if K_p<1:
1568           :param safety_factor: safety factor for addjustment of pressure and velocity tolerance from stopping criteria                 raise ValueError("K_p need to be greater or equal to 1.")
          :type safety_factor: ``float``  
          """  
          if not gamma == None:  
             if gamma<=0 or gamma>=1:  
                raise ValueError,"Initial gamma needs to be positive and less than 1."  
          else:  
             gamma=self.__gamma  
   
          if not gamma_min == None:  
             if gamma_min<=0 or gamma_min>=1:  
                raise ValueError,"gamma_min needs to be positive and less than 1."  
          else:  
             gamma_min = self.__gamma_min  
   
          if not chi_max == None:  
             if chi_max<=0 or chi_max>=1:  
                raise ValueError,"chi_max needs to be positive and less than 1."  
          else:  
             chi_max = self.__chi_max  
   
          if not omega_div == None:  
             if omega_div<=0 or omega_div >=1:  
                raise ValueError,"omega_div needs to be positive and less than 1."  
1569           else:           else:
1570              omega_div=self.__omega_div              K_p=self.__K_p
1571    
1572           if not omega_conv == None:           if not K_v == None:
1573              if omega_conv<1:              if K_v<1:
1574                 raise ValueError,"omega_conv needs to be greater or equal to 1."                 raise ValueError("K_v need to be greater or equal to 1.")
1575           else:           else:
1576              omega_conv=self.__omega_conv              K_v=self.__K_v
   
          if not rtol_min == None:  
             if rtol_min<=0 or rtol_min>=1:  
                raise ValueError,"rtol_min needs to be positive and less than 1."  
          else:  
             rtol_min=self.__rtol_min  
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 chi == None:           if not rtol_min == None:
1585              if chi<=0 or chi>1:              if rtol_min<=0 or rtol_min>=1:
1586                 raise ValueError,"chi needs to be positive and less than 1."                 raise ValueError("rtol_min needs to be positive and less than 1.")
1587           else:           else:
1588              chi=self.__chi              rtol_min=self.__rtol_min
1589    
1590           if not C_p == None:           if not chi_max == None:
1591              if C_p<1:              if chi_max<=0 or chi_max>=1:
1592                 raise ValueError,"C_p need to be greater or equal to 1."                 raise ValueError("chi_max needs to be positive and less than 1.")
1593           else:           else:
1594              C_p=self.__C_p              chi_max = self.__chi_max
1595    
1596           if not C_v == None:           if not reduction_factor == None:
1597              if C_v<1:              if reduction_factor<=0 or reduction_factor>1:
1598                 raise ValueError,"C_v need to be greater or equal to 1."                 raise ValueError("reduction_factor need to be between zero and one.")
1599           else:           else:
1600              C_v=self.__C_v              reduction_factor=self.__reduction_factor
1601    
1602           if not safety_factor == None:           if not theta == None:
1603              if safety_factor<=0 or safety_factor>1:              if theta<=0 or theta>1:
1604                 raise ValueError,"safety_factor need to be between zero and one."                 raise ValueError("theta need to be between zero and one.")
1605           else:           else:
1606              safety_factor=self.__safety_factor              theta=self.__theta
1607    
1608           if gamma<gamma_min:           if rtol_min>=rtol_max:
1609                 raise ValueError,"gamma = %e needs to be greater or equal gamma_min = %e."%(gamma,gamma_min)               raise ValueError("rtol_max = %e needs to be greater than rtol_min = %e"%(rtol_max,rtol_min))
          if rtol_max<=rtol_min:  
                raise ValueError,"rtol_max = %e needs to be greater rtol_min = %e."%(rtol_max,rtol_min)  
               
          self.__gamma = gamma  
          self.__gamma_min = gamma_min  
1610           self.__chi_max = chi_max           self.__chi_max = chi_max
          self.__omega_div = omega_div  
          self.__omega_conv = omega_conv  
          self.__rtol_min = rtol_min  
1611           self.__rtol_max = rtol_max           self.__rtol_max = rtol_max
1612           self.__chi = chi           self.__K_p = K_p
1613           self.__C_p = C_p           self.__K_v = K_v
1614           self.__C_v = C_v           self.__reduction_factor = reduction_factor
1615           self.__safety_factor = safety_factor           self.__theta = theta
1616             self.__rtol_min=rtol_min
1617    
1618        #=============================================================        #=============================================================
1619        def inner_pBv(self,p,Bv):        def inner_pBv(self,p,Bv):
# Line 1657  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 1668  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 1678  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 1688  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 1698  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 1707  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 1717  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 1726  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 1753  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 1787  class HomogeneousSaddlePointProblem(obje Line 1756  class HomogeneousSaddlePointProblem(obje
1756           self.verbose=verbose           self.verbose=verbose
1757           rtol=self.getTolerance()           rtol=self.getTolerance()
1758           atol=self.getAbsoluteTolerance()           atol=self.getAbsoluteTolerance()
1759    
1760             K_p=self.__K_p
1761             K_v=self.__K_v
1762           correction_step=0           correction_step=0
1763           converged=False           converged=False
          error=None  
1764           chi=None           chi=None
1765           gamma=self.__gamma           eps=None
1766           C_p=self.__C_p  
1767           C_v=self.__C_v           if self.verbose: print(("HomogeneousSaddlePointProblem: start iteration: rtol= %e, atol=%e"%(rtol, atol)))
1768           while not converged:           while not converged:
1769                if error== None or chi == None:  
1770                    gamma_new=gamma/self.__omega_conv               # get tolerance for velecity increment:
1771                else:               if chi == None:
1772                   if chi < self.__chi_max:                  rtol_v=self.__rtol_max
1773                      gamma_new=min(max(gamma*self.__omega_conv,1-chi*error/(self.__safety_factor*ATOL)), 1-chi/self.__chi_max)               else:
1774                   else:                  rtol_v=min(chi/K_v,self.__rtol_max)
1775                      gamma_new=gamma*self.__omega_div               rtol_v=max(rtol_v, self.__rtol_min)
1776                gamma=max(gamma_new, self.__gamma_min)               if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: rtol_v= %e"%(correction_step,rtol_v)))
1777                # calculate velocity for current pressure:               # get velocity increment:
1778                rtol_v=min(max(gamma/(1.+gamma)/C_v,self.__rtol_min), self.__rtol_max)               dv1=self.getDV(p,v,rtol_v)
1779                rtol_p=min(max(gamma/C_p, self.__rtol_min), self.__rtol_max)               v1=v+dv1
1780                self.__subtol=rtol_p**2               Bv1=self.Bv(v1, rtol_v)
1781                if self.verbose: print "HomogeneousSaddlePointProblem: step %s: gamma = %e, rtol_v= %e, rtol_p=%e"%(correction_step,gamma,rtol_v,rtol_p)               norm_Bv1=self.norm_Bv(Bv1)
1782                if self.verbose: print "HomogeneousSaddlePointProblem: subtolerance: %e"%self.__subtol               norm_dv1=self.norm_v(dv1)
1783                # calculate velocity for current pressure: A*dv= F-A*v-B*p               if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: norm_Bv1 = %e, norm_dv1 = %e"%(correction_step, norm_Bv1, norm_dv1)))
1784                dv1=self.getDV(p,v,rtol_v)               if norm_dv1*self.__theta < norm_Bv1:
1785                v1=v+dv1                  # get tolerance for pressure increment:
1786                Bv1=self.Bv(v1, self.__subtol)                  large_Bv1=True
1787                norm_Bv1=self.norm_Bv(Bv1)                  if chi == None or eps == None:
1788                norm_dv1=self.norm_v(dv1)                     rtol_p=self.__rtol_max
1789                norm_v1=self.norm_v(v1)                  else:
1790                ATOL=norm_v1*rtol+atol                     rtol_p=min(chi**2*eps/K_p/norm_Bv1, self.__rtol_max)
1791                if self.verbose: print "HomogeneousSaddlePointProblem: step %s: Bv = %e, dv = %e, v=%e"%(correction_step,norm_Bv1, norm_dv1, norm_v1)                  self.__subtol=max(rtol_p**2, self.__rtol_min)
1792                if not ATOL>0: raise ValueError,"overall absolute tolerance needs to be positive."                  if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: rtol_p= %e"%(correction_step,rtol_p)))
1793                if max(norm_Bv1, norm_dv1) <= ATOL:                  # now we solve for the pressure increment dp from B*A^{-1}B^* dp = Bv1
1794                    converged=True                  if usePCG:
                   v=v1  
               else:  
                   # now we solve for the pressure increment dp from B*A^{-1}B^* dp = Bv1  
                   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)
1796                      v2=r[0]                      v2=r[0]
1797                      Bv2=r[1]                      Bv2=r[1]
1798                    else:                  else:
1799                        # don't use!!!!
1800                      dp=GMRES(self.solve_prec(Bv1,self.__subtol),self.__Aprod_GMRES, 0*p, self.__inner_GMRES,atol=0, rtol=rtol_p,iter_max=max_iter, iter_restart=iter_restart, verbose=self.verbose)                      dp=GMRES(self.solve_prec(Bv1,self.__subtol),self.__Aprod_GMRES, 0*p, self.__inner_GMRES,atol=0, rtol=rtol_p,iter_max=max_iter, iter_restart=iter_restart, verbose=self.verbose)
1801                      dv2=self.solve_AinvBt(dp, self.__subtol)                      dv2=self.solve_AinvBt(dp, self.__subtol)
1802                      v2=v1-dv2                      v2=v1-dv2
1803                      Bv2=self.Bv(v2, self.__subtol)                      Bv2=self.Bv(v2, self.__subtol)
1804                    #                  p2=p+dp
1805                    # convergence indicators:               else:
1806                    #                  large_Bv1=False
1807                    norm_v2=self.norm_v(v2)                  v2=v1
1808                    norm_dv2=self.norm_v(v2-v)                  p2=p
1809                    error_new=max(norm_dv2, norm_Bv1)               # update business:
1810                    norm_Bv2=self.norm_Bv(Bv2)               norm_dv2=self.norm_v(v2-v)
1811                    if self.verbose: print "HomogeneousSaddlePointProblem: step %s (part 2): Bv = %e, dv = %e, v=%e"%(correction_step,norm_Bv2, norm_dv2, norm_v2)               norm_v2=self.norm_v(v2)
1812                    if error !=None:               if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: v2 = %e, norm_dv2 = %e"%(correction_step, norm_v2, self.norm_v(v2-v))))
1813                        chi_new=error_new/error               eps, eps_old = max(norm_Bv1, norm_dv2), eps
1814                        if self.verbose: print "HomogeneousSaddlePointProblem: step %s: convergence rate = %e, est. error = %e"%(correction_step,chi_new, error_new)               if eps_old == None:
1815                        if chi != None:                    chi, chi_old = None, chi
1816                            gamma0=max(gamma, 1-chi/chi_new)               else:
1817                            C_p*=gamma0/gamma                    chi, chi_old = min(eps/ eps_old, self.__chi_max), chi
1818                            C_v*=gamma0/gamma*(1+gamma)/(1+gamma0)               if eps != None:
1819                        chi=chi_new                   if chi !=None:
1820                    else:                      if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: convergence rate = %e, correction = %e"%(correction_step,chi, eps)))
1821                        if self.verbose: print "HomogeneousSaddlePointProblem: step %s: est. error = %e"%(correction_step, error_new)                   else:
1822                        if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: correction = %e"%(correction_step, eps)))
1823                    error = error_new               if eps <= rtol*norm_v2+atol :
1824                    correction_step+=1                   converged = True
1825                    if correction_step>max_correction_steps:               else:
1826                          raise CorrectionFailed,"Given up after %d correction steps."%correction_step                   if correction_step>=max_correction_steps:
1827                    v,p=v2,p+dp                        raise CorrectionFailed("Given up after %d correction steps."%correction_step)
1828           if self.verbose: print "HomogeneousSaddlePointProblem: tolerance reached after %s steps."%correction_step                   if chi_old!=None:
1829       return v,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)
1831                        if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: new adjustment factor K = %e"%(correction_step,K_p)))
1832                 correction_step+=1
1833                 v,p =v2, p2
1834             if self.verbose: print(("HomogeneousSaddlePointProblem: tolerance reached after %s steps."%correction_step))
1835             return v,p
1836        #========================================================================        #========================================================================
1837        def setTolerance(self,tolerance=1.e-4):        def setTolerance(self,tolerance=1.e-4):
1838           """           """
# Line 1869  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 1889  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.2793  
changed lines
  Added in v.3911

  ViewVC Help
Powered by ViewVC 1.1.26