/[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 3981 by jfenwick, Fri Sep 21 02:47:54 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)  # http://www.uq.edu.au
 # http://www.uq.edu.au/esscc  
6  #  #
7  # Primary Business: Queensland, Australia  # Primary Business: Queensland, Australia
8  # Licensed under the Open Software License version 3.0  # Licensed under the Open Software License version 3.0
9  # http://www.opensource.org/licenses/osl-3.0.php  # http://www.opensource.org/licenses/osl-3.0.php
10  #  #
11  ########################################################  # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12    # Development since 2012 by School of Earth Sciences
13    #
14    ##############################################################################
15    
16  __copyright__="""Copyright (c) 2003-2009 by University of Queensland  __copyright__="""Copyright (c) 2003-2012 by University of Queensland
17  Earth Systems Science Computational Center (ESSCC)  http://www.uq.edu.au
 http://www.uq.edu.au/esscc  
18  Primary Business: Queensland, Australia"""  Primary Business: Queensland, Australia"""
19  __license__="""Licensed under the Open Software License version 3.0  __license__="""Licensed under the Open Software License version 3.0
20  http://www.opensource.org/licenses/osl-3.0.php"""  http://www.opensource.org/licenses/osl-3.0.php"""
# Line 39  Currently includes: Line 40  Currently includes:
40  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
41    
42    
43  import escript  from . import escript
44  import linearPDEs  from . import linearPDEs
45  import numpy  import numpy
46  import util  from . import util
47  import math  import math
48    
49  class TimeIntegrationManager:  class TimeIntegrationManager:
# Line 66  class TimeIntegrationManager: Line 67  class TimeIntegrationManager:
67       Sets up the value manager where ``inital_values`` are the initial values       Sets up the value manager where ``inital_values`` are the initial values
68       and p is the order used for extrapolation.       and p is the order used for extrapolation.
69       """       """
70       if kwargs.has_key("p"):       if "p" in kwargs:
71              self.__p=kwargs["p"]              self.__p=kwargs["p"]
72       else:       else:
73              self.__p=1              self.__p=1
74       if kwargs.has_key("time"):       if "time" in kwargs:
75              self.__t=kwargs["time"]              self.__t=kwargs["time"]
76       else:       else:
77              self.__t=0.              self.__t=0.
# Line 153  class Projector: Line 154  class Projector:
154      """      """
155      return self.__pde.getSolverOptions()      return self.__pde.getSolverOptions()
156    
157      def getValue(self, input_data):
158        """
159        Projects ``input_data`` onto a continuous function.
160    
161        :param input_data: the data to be projected
162        """
163        return self(input_data)
164    
165    def __call__(self, input_data):    def __call__(self, input_data):
166      """      """
167      Projects ``input_data`` onto a continuous function.      Projects ``input_data`` onto a continuous function.
# Line 291  class NoPDE: Line 300  class NoPDE:
300           """           """
301           if self.__u==None:           if self.__u==None:
302              if self.__D==None:              if self.__D==None:
303                 raise ValueError,"coefficient D is undefined"                 raise ValueError("coefficient D is undefined")
304              D=escript.Data(self.__D,self.__function_space)              D=escript.Data(self.__D,self.__function_space)
305              if D.getRank()>1:              if D.getRank()>1:
306                 raise ValueError,"coefficient D must have rank 0 or 1"                 raise ValueError("coefficient D must have rank 0 or 1")
307              if self.__Y==None:              if self.__Y==None:
308                 self.__u=escript.Data(0.,D.getShape(),self.__function_space)                 self.__u=escript.Data(0.,D.getShape(),self.__function_space)
309              else:              else:
310                 self.__u=util.quotient(self.__Y,D)                 self.__u=1./D*self.__Y
311              if not self.__q==None:              if not self.__q==None:
312                  q=util.wherePositive(escript.Data(self.__q,self.__function_space))                  q=util.wherePositive(escript.Data(self.__q,self.__function_space))
313                  self.__u*=(1.-q)                  self.__u*=(1.-q)
# Line 332  class Locator: Line 341  class Locator:
341         iterative=False         iterative=False
342         if isinstance(x, list):         if isinstance(x, list):
343             if len(x)==0:             if len(x)==0:
344                raise "ValueError", "At least one point must be given."                raise ValueError("At least one point must be given.")
345             try:             try:
346               iter(x[0])               iter(x[0])
347               iterative=True               iterative=True
348             except TypeError:             except TypeError:
349               iterative=False               iterative=False
350           xxx=self.__function_space.getX()
351         if iterative:         if iterative:
352             self.__id=[]             self.__id=[]
353             for p in x:             for p in x:
354                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())
355         else:         else:
356             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()
357    
358       def __str__(self):       def __str__(self):
359         """         """
# Line 421  class Locator: Line 431  class Locator:
431                  return out                  return out
432          else:          else:
433             return data             return data
434              
435         def setValue(self, data, v):
436          """
437          Sets the value of the ``data`` at the Locator.
438          """
439          if isinstance(data, escript.Data):
440             if data.getFunctionSpace()!=self.getFunctionSpace():
441               raise TypeError, "setValue: FunctionSpace of Locator and Data object must match."
442             data.expand()  
443             id=self.getId()
444             if isinstance(id, list):
445              for i in id:
446               data._setTupleForGlobalDataPoint(i[1], i[0], v)
447             else:
448               data._setTupleForGlobalDataPoint(id[1], id[0], v)
449          else:
450               raise TypeError, "setValue: Invalid argument type."
451    
452    
453  def getInfLocator(arg):  def getInfLocator(arg):
# Line 428  def getInfLocator(arg): Line 455  def getInfLocator(arg):
455      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.
456      """      """
457      if not isinstance(arg, escript.Data):      if not isinstance(arg, escript.Data):
458      raise TypeError,"getInfLocator: Unknown argument type."         raise TypeError("getInfLocator: Unknown argument type.")
459      a_inf=util.inf(arg)      a_inf=util.inf(arg)
460      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
461      x=arg.getFunctionSpace().getX()      x=arg.getFunctionSpace().getX()
# Line 440  def getSupLocator(arg): Line 467  def getSupLocator(arg):
467      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.
468      """      """
469      if not isinstance(arg, escript.Data):      if not isinstance(arg, escript.Data):
470      raise TypeError,"getInfLocator: Unknown argument type."         raise TypeError("getInfLocator: Unknown argument type.")
471      a_inf=util.sup(arg)      a_inf=util.sup(arg)
472      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
473      x=arg.getFunctionSpace().getX()      x=arg.getFunctionSpace().getX()
# Line 539  def PCG(r, Aprod, x, Msolve, bilinearfor Line 566  def PCG(r, Aprod, x, Msolve, bilinearfor
566     rhat=Msolve(r)     rhat=Msolve(r)
567     d = rhat     d = rhat
568     rhat_dot_r = bilinearform(rhat, r)     rhat_dot_r = bilinearform(rhat, r)
569     if rhat_dot_r<0: raise NegativeNorm,"negative norm."     if rhat_dot_r<0: raise NegativeNorm("negative norm.")
570     norm_r0=math.sqrt(rhat_dot_r)     norm_r0=math.sqrt(rhat_dot_r)
571     atol2=atol+rtol*norm_r0     atol2=atol+rtol*norm_r0
572     if atol2<=0:     if atol2<=0:
573        raise ValueError,"Non-positive tolarance."        raise ValueError("Non-positive tolarance.")
574     atol2=max(atol2, 100. * util.EPSILON * norm_r0)     atol2=max(atol2, 100. * util.EPSILON * norm_r0)
575    
576     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)))
577    
578    
579     while not math.sqrt(rhat_dot_r) <= atol2:     while not math.sqrt(rhat_dot_r) <= atol2:
580         iter+=1         iter+=1
581         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)
582    
583         q=Aprod(d)         q=Aprod(d)
584         alpha = rhat_dot_r / bilinearform(d, q)         alpha = rhat_dot_r / bilinearform(d, q)
585         x += alpha * d         x += alpha * d
586         if isinstance(q,ArithmeticTuple):         if isinstance(q,ArithmeticTuple):
587         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__
588         else:         else:
589             r += (-alpha) * q             r += (-alpha) * q
590         rhat=Msolve(r)         rhat=Msolve(r)
# Line 567  def PCG(r, Aprod, x, Msolve, bilinearfor Line 594  def PCG(r, Aprod, x, Msolve, bilinearfor
594         d=rhat         d=rhat
595    
596         rhat_dot_r = rhat_dot_r_new         rhat_dot_r = rhat_dot_r_new
597         if rhat_dot_r<0: raise NegativeNorm,"negative norm."         if rhat_dot_r<0: raise NegativeNorm("negative norm.")
598         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))))
599     if verbose: print "PCG: tolerance reached after %s steps."%iter     if verbose: print(("PCG: tolerance reached after %s steps."%iter))
600     return x,r,math.sqrt(rhat_dot_r)     return x,r,math.sqrt(rhat_dot_r)
601    
602  class Defect(object):  class Defect(object):
# Line 603  class Defect(object): Line 630  class Defect(object):
630          :note: by default ``sqrt(self.bilinearform(x,x)`` is returned.          :note: by default ``sqrt(self.bilinearform(x,x)`` is returned.
631          """          """
632          s=self.bilinearform(x,x)          s=self.bilinearform(x,x)
633          if s<0: raise NegativeNorm,"negative norm."          if s<0: raise NegativeNorm("negative norm.")
634          return math.sqrt(s)          return math.sqrt(s)
635    
636      def eval(self,x):      def eval(self,x):
# Line 627  class Defect(object): Line 654  class Defect(object):
654          :param inc: relative increment length          :param inc: relative increment length
655          :type inc: positive ``float``          :type inc: positive ``float``
656          """          """
657          if inc<=0: raise ValueError,"positive increment required."          if inc<=0: raise ValueError("positive increment required.")
658          self.__inc=inc          self.__inc=inc
659    
660      def getDerivativeIncrementLength(self):      def getDerivativeIncrementLength(self):
# Line 698  def NewtonGMRES(defect, x, iter_max=100, Line 725  def NewtonGMRES(defect, x, iter_max=100,
725     :rtype: same type as the initial guess     :rtype: same type as the initial guess
726     """     """
727     lmaxit=iter_max     lmaxit=iter_max
728     if atol<0: raise ValueError,"atol needs to be non-negative."     if atol<0: raise ValueError("atol needs to be non-negative.")
729     if rtol<0: raise ValueError,"rtol needs to be non-negative."     if rtol<0: raise ValueError("rtol needs to be non-negative.")
730     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.")
731     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)
732     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)
733    
734     F=defect(x)     F=defect(x)
735     fnrm=defect.norm(F)     fnrm=defect.norm(F)
736     stop_tol=atol + rtol*fnrm     stop_tol=atol + rtol*fnrm
737     subtol=subtol_max     subtol=subtol_max
738     if verbose: print "NewtonGMRES: initial residual = %e."%fnrm     if verbose: print(("NewtonGMRES: initial residual = %e."%fnrm))
739     if verbose: print "             tolerance = %e."%subtol     if verbose: print(("             tolerance = %e."%subtol))
740     iter=1     iter=1
741     #     #
742     # main iteration loop     # main iteration loop
743     #     #
744     while not fnrm<=stop_tol:     while not fnrm<=stop_tol:
745              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)
746              #              #
747          #   adjust subtol_          #   adjust subtol_
748          #          #
749              if iter > 1:              if iter > 1:
750             rat=fnrm/fnrmo                 rat=fnrm/fnrmo
751                 subtol_old=subtol                 subtol_old=subtol
752             subtol=gamma*rat**2                 subtol=gamma*rat**2
753             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)
754             subtol=max(min(subtol,subtol_max), .5*stop_tol/fnrm)                 subtol=max(min(subtol,subtol_max), .5*stop_tol/fnrm)
755          #          #
756          # calculate newton increment xc          # calculate newton increment xc
757              #     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 759  def NewtonGMRES(defect, x, iter_max=100,
759              #     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
760              #              #
761              #              #
762              if verbose: print "             subiteration (GMRES) is called with relative tolerance %e."%subtol              if verbose: print(("             subiteration (GMRES) is called with relative tolerance %e."%subtol))
763              try:              try:
764                 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)
765              except MaxIterReached:              except MaxIterReached:
766                 raise MaxIterReached,"maximum number of %s steps reached."%iter_max                 raise MaxIterReached("maximum number of %s steps reached."%iter_max)
767              if sub_iter<0:              if sub_iter<0:
768                 iter+=sub_iter_max                 iter+=sub_iter_max
769              else:              else:
770                 iter+=sub_iter                 iter+=sub_iter
771              # ====              # ====
772          x+=xc              x+=xc
773              F=defect(x)              F=defect(x)
774          iter+=1              iter+=1
775              fnrmo, fnrm=fnrm, defect.norm(F)              fnrmo, fnrm=fnrm, defect.norm(F)
776              if verbose: print "             step %s: residual %e."%(iter,fnrm)              if verbose: print(("             step %s: residual %e."%(iter,fnrm)))
777     if verbose: print "NewtonGMRES: completed after %s steps."%iter     if verbose: print(("NewtonGMRES: completed after %s steps."%iter))
778     return x     return x
779    
780  def __givapp(c,s,vin):  def __givapp(c,s,vin):
# Line 763  def __givapp(c,s,vin): Line 790  def __givapp(c,s,vin):
790      else:      else:
791          for i in range(len(c)):          for i in range(len(c)):
792              w1=c[i]*vrot[i]-s[i]*vrot[i+1]              w1=c[i]*vrot[i]-s[i]*vrot[i+1]
793          w2=s[i]*vrot[i]+c[i]*vrot[i+1]              w2=s[i]*vrot[i]+c[i]*vrot[i+1]
794              vrot[i]=w1              vrot[i]=w1
795              vrot[i+1]=w2              vrot[i+1]=w2
796      return vrot      return vrot
# Line 783  def __FDGMRES(F0, defect, x0, atol, iter Line 810  def __FDGMRES(F0, defect, x0, atol, iter
810     iter=0     iter=0
811     while rho > atol and iter<iter_restart-1:     while rho > atol and iter<iter_restart-1:
812          if iter  >= iter_max:          if iter  >= iter_max:
813              raise MaxIterReached,"maximum number of %s steps reached."%iter_max              raise MaxIterReached("maximum number of %s steps reached."%iter_max)
814    
815          p=defect.derivative(F0,x0,v[iter], v_is_normalised=True)          p=defect.derivative(F0,x0,v[iter], v_is_normalised=True)
816          v.append(p)          v.append(p)
# Line 845  def __FDGMRES(F0, defect, x0, atol, iter Line 872  def __FDGMRES(F0, defect, x0, atol, iter
872            i=i-1            i=i-1
873       xhat=v[iter-1]*y[iter-1]       xhat=v[iter-1]*y[iter-1]
874       for i in range(iter-1):       for i in range(iter-1):
875      xhat += v[i]*y[i]         xhat += v[i]*y[i]
876     else :     else :
877        xhat=v[0] * 0        xhat=v[0] * 0
878    
# Line 904  def GMRES(r, Aprod, x, bilinearform, ato Line 931  def GMRES(r, Aprod, x, bilinearform, ato
931     iter=0     iter=0
932     if rtol>0:     if rtol>0:
933        r_dot_r = bilinearform(r, r)        r_dot_r = bilinearform(r, r)
934        if r_dot_r<0: raise NegativeNorm,"negative norm."        if r_dot_r<0: raise NegativeNorm("negative norm.")
935        atol2=atol+rtol*math.sqrt(r_dot_r)        atol2=atol+rtol*math.sqrt(r_dot_r)
936        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)))
937     else:     else:
938        atol2=atol        atol2=atol
939        if verbose: print "GMRES: absolute tolerance = %e"%atol2        if verbose: print(("GMRES: absolute tolerance = %e"%atol2))
940     if atol2<=0:     if atol2<=0:
941        raise ValueError,"Non-positive tolarance."        raise ValueError("Non-positive tolarance.")
942    
943     while True:     while True:
944        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)
945        if restarted:        if restarted:
946           r2 = r-Aprod(x-x2)           r2 = r-Aprod(x-x2)
947        else:        else:
# Line 923  def GMRES(r, Aprod, x, bilinearform, ato Line 950  def GMRES(r, Aprod, x, bilinearform, ato
950        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)
951        iter+=iter_restart        iter+=iter_restart
952        if stopped: break        if stopped: break
953        if verbose: print "GMRES: restart."        if verbose: print("GMRES: restart.")
954        restarted=True        restarted=True
955     if verbose: print "GMRES: tolerance has been reached."     if verbose: print("GMRES: tolerance has been reached.")
956     return x     return x
957    
958  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 965  def _GMRESm(r, Aprod, x, bilinearform, a
965     v=[]     v=[]
966    
967     r_dot_r = bilinearform(r, r)     r_dot_r = bilinearform(r, r)
968     if r_dot_r<0: raise NegativeNorm,"negative norm."     if r_dot_r<0: raise NegativeNorm("negative norm.")
969     rho=math.sqrt(r_dot_r)     rho=math.sqrt(r_dot_r)
970    
971     v.append(r/rho)     v.append(r/rho)
972     g[0]=rho     g[0]=rho
973    
974     if verbose: print "GMRES: initial residual %e (absolute tolerance = %e)"%(rho,atol)     if verbose: print(("GMRES: initial residual %e (absolute tolerance = %e)"%(rho,atol)))
975     while not (rho<=atol or iter==iter_restart):     while not (rho<=atol or iter==iter_restart):
976    
977      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)
978    
979          if P_R!=None:          if P_R!=None:
980              p=Aprod(P_R(v[iter]))              p=Aprod(P_R(v[iter]))
981          else:          else:
982          p=Aprod(v[iter])              p=Aprod(v[iter])
983      v.append(p)          v.append(p)
984    
985      v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))          v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
986    
987  # Modified Gram-Schmidt  # Modified Gram-Schmidt
988      for j in range(iter+1):          for j in range(iter+1):
989        h[j,iter]=bilinearform(v[j],v[iter+1])            h[j,iter]=bilinearform(v[j],v[iter+1])
990        v[iter+1]-=h[j,iter]*v[j]            v[iter+1]-=h[j,iter]*v[j]
991    
992      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]))
993      v_norm2=h[iter+1,iter]          v_norm2=h[iter+1,iter]
994    
995  # Reorthogonalize if needed  # Reorthogonalize if needed
996      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)
997       for j in range(iter+1):           for j in range(iter+1):
998          hr=bilinearform(v[j],v[iter+1])              hr=bilinearform(v[j],v[iter+1])
999              h[j,iter]=h[j,iter]+hr              h[j,iter]=h[j,iter]+hr
1000              v[iter+1] -= hr*v[j]              v[iter+1] -= hr*v[j]
1001    
1002       v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))           v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
1003       h[iter+1,iter]=v_norm2           h[iter+1,iter]=v_norm2
1004    
1005  #   watch out for happy breakdown  #   watch out for happy breakdown
1006          if not v_norm2 == 0:          if not v_norm2 == 0:
1007           v[iter+1]=v[iter+1]/h[iter+1,iter]           v[iter+1]=v[iter+1]/h[iter+1,iter]
1008    
1009  #   Form and store the information for the new Givens rotation  #   Form and store the information for the new Givens rotation
1010      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])
1011      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])
1012    
1013      if mu!=0 :          if mu!=0 :
1014          c[iter]=h[iter,iter]/mu                  c[iter]=h[iter,iter]/mu
1015          s[iter]=-h[iter+1,iter]/mu                  s[iter]=-h[iter+1,iter]/mu
1016          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]
1017          h[iter+1,iter]=0.0                  h[iter+1,iter]=0.0
1018                  gg=__givapp(c[iter],s[iter],[g[iter],g[iter+1]])                  gg=__givapp(c[iter],s[iter],[g[iter],g[iter+1]])
1019                  g[iter]=gg[0]                  g[iter]=gg[0]
1020                  g[iter+1]=gg[1]                  g[iter+1]=gg[1]
1021  # Update the residual norm  # Update the residual norm
1022    
1023          rho=abs(g[iter+1])          rho=abs(g[iter+1])
1024          if verbose: print "GMRES: iteration step %s: residual %e"%(iter,rho)          if verbose: print(("GMRES: iteration step %s: residual %e"%(iter,rho)))
1025      iter+=1          iter+=1
1026    
1027  # At this point either iter > iter_max or rho < tol.  # At this point either iter > iter_max or rho < tol.
1028  # It's time to compute x and leave.  # It's time to compute x and leave.
1029    
1030     if verbose: print "GMRES: iteration stopped after %s step."%iter     if verbose: print(("GMRES: iteration stopped after %s step."%iter))
1031     if iter > 0 :     if iter > 0 :
1032       y=numpy.zeros(iter,numpy.float64)       y=numpy.zeros(iter,numpy.float64)
1033       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 1038  def _GMRESm(r, Aprod, x, bilinearform, a
1038            i=i-1            i=i-1
1039       xhat=v[iter-1]*y[iter-1]       xhat=v[iter-1]*y[iter-1]
1040       for i in range(iter-1):       for i in range(iter-1):
1041      xhat += v[i]*y[i]         xhat += v[i]*y[i]
1042     else:     else:
1043       xhat=v[0] * 0       xhat=v[0] * 0
1044     if P_R!=None:     if P_R!=None:
# Line 1084  def MINRES(r, Aprod, x, Msolve, bilinear Line 1111  def MINRES(r, Aprod, x, Msolve, bilinear
1111      y = Msolve(r)      y = Msolve(r)
1112      beta1 = bilinearform(y,r)      beta1 = bilinearform(y,r)
1113    
1114      if beta1< 0: raise NegativeNorm,"negative norm."      if beta1< 0: raise NegativeNorm("negative norm.")
1115    
1116      #  If r = 0 exactly, stop with x      #  If r = 0 exactly, stop with x
1117      if beta1==0: return x      if beta1==0: return x
# Line 1119  def MINRES(r, Aprod, x, Msolve, bilinear Line 1146  def MINRES(r, Aprod, x, Msolve, bilinear
1146      # --------------------------------------------------------------------      # --------------------------------------------------------------------
1147      while not rnorm<=atol+rtol*Anorm*ynorm:    #  checks ||r|| < (||A|| ||x||) * TOL      while not rnorm<=atol+rtol*Anorm*ynorm:    #  checks ||r|| < (||A|| ||x||) * TOL
1148    
1149      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)
1150          iter    = iter  +  1          iter    = iter  +  1
1151    
1152          #-----------------------------------------------------------------          #-----------------------------------------------------------------
# Line 1149  def MINRES(r, Aprod, x, Msolve, bilinear Line 1176  def MINRES(r, Aprod, x, Msolve, bilinear
1176          y = Msolve(r2)          y = Msolve(r2)
1177          oldb   = beta                           # oldb = betak          oldb   = beta                           # oldb = betak
1178          beta   = bilinearform(y,r2)             # beta = betak+1^2          beta   = bilinearform(y,r2)             # beta = betak+1^2
1179          if beta < 0: raise NegativeNorm,"negative norm."          if beta < 0: raise NegativeNorm("negative norm.")
1180    
1181          beta   = math.sqrt( beta )          beta   = math.sqrt( beta )
1182          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 1283  def TFQMR(r, Aprod, x, bilinearform, ato
1283    theta = 0.0;    theta = 0.0;
1284    eta = 0.0;    eta = 0.0;
1285    rho=bilinearform(r,r)    rho=bilinearform(r,r)
1286    if rho < 0: raise NegativeNorm,"negative norm."    if rho < 0: raise NegativeNorm("negative norm.")
1287    tau = math.sqrt(rho)    tau = math.sqrt(rho)
1288    norm_r0=tau    norm_r0=tau
1289    while tau>atol+rtol*norm_r0:    while tau>atol+rtol*norm_r0:
1290      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)
1291    
1292      sigma = bilinearform(r,v)      sigma = bilinearform(r,v)
1293      if sigma == 0.0: raise IterationBreakDown,'TFQMR breakdown, sigma=0'      if sigma == 0.0: raise IterationBreakDown('TFQMR breakdown, sigma=0')
1294    
1295      alpha = rho / sigma      alpha = rho / sigma
1296    
# Line 1291  def TFQMR(r, Aprod, x, bilinearform, ato Line 1318  def TFQMR(r, Aprod, x, bilinearform, ato
1318  #  #
1319  #  Try to terminate the iteration at each pass through the loop  #  Try to terminate the iteration at each pass through the loop
1320  #  #
1321      if rho == 0.0: raise IterationBreakDown,'TFQMR breakdown, rho=0'      if rho == 0.0: raise IterationBreakDown('TFQMR breakdown, rho=0')
1322    
1323      rhon = bilinearform(r,w)      rhon = bilinearform(r,w)
1324      beta = rhon / rho;      beta = rhon / rho;
# Line 1363  class ArithmeticTuple(object): Line 1390  class ArithmeticTuple(object):
1390         try:         try:
1391             l=len(other)             l=len(other)
1392             if l!=len(self):             if l!=len(self):
1393                 raise ValueError,"length of arguments don't match."                 raise ValueError("length of arguments don't match.")
1394             for i in range(l): out.append(self[i]*other[i])             for i in range(l): out.append(self[i]*other[i])
1395         except TypeError:         except TypeError:
1396             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 1409  class ArithmeticTuple(object):
1409         try:         try:
1410             l=len(other)             l=len(other)
1411             if l!=len(self):             if l!=len(self):
1412                 raise ValueError,"length of arguments don't match."                 raise ValueError("length of arguments don't match.")
1413             for i in range(l): out.append(other[i]*self[i])             for i in range(l): out.append(other[i]*self[i])
1414         except TypeError:         except TypeError:
1415             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 1439  class ArithmeticTuple(object):
1439         try:         try:
1440             l=len(other)             l=len(other)
1441             if l!=len(self):             if l!=len(self):
1442                 raise ValueError,"length of arguments don't match."                 raise ValueError("length of arguments don't match.")
1443             for i in range(l): out.append(other[i]/self[i])             for i in range(l): out.append(other[i]/self[i])
1444         except TypeError:         except TypeError:
1445             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 1453  class ArithmeticTuple(object):
1453         :type other: ``ArithmeticTuple``         :type other: ``ArithmeticTuple``
1454         """         """
1455         if len(self) != len(other):         if len(self) != len(other):
1456             raise ValueError,"tuple lengths must match."             raise ValueError("tuple lengths must match.")
1457         for i in range(len(self)):         for i in range(len(self)):
1458             self.__items[i]+=other[i]             self.__items[i]+=other[i]
1459         return self         return self
# Line 1442  class ArithmeticTuple(object): Line 1469  class ArithmeticTuple(object):
1469         try:         try:
1470             l=len(other)             l=len(other)
1471             if l!=len(self):             if l!=len(self):
1472                 raise ValueError,"length of arguments don't match."                 raise ValueError("length of arguments don't match.")
1473             for i in range(l): out.append(self[i]+other[i])             for i in range(l): out.append(self[i]+other[i])
1474         except TypeError:         except TypeError:
1475             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 1486  class ArithmeticTuple(object):
1486         try:         try:
1487             l=len(other)             l=len(other)
1488             if l!=len(self):             if l!=len(self):
1489                 raise ValueError,"length of arguments don't match."                 raise ValueError("length of arguments don't match.")
1490             for i in range(l): out.append(self[i]-other[i])             for i in range(l): out.append(self[i]-other[i])
1491         except TypeError:         except TypeError:
1492             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 1500  class ArithmeticTuple(object):
1500         :type other: ``ArithmeticTuple``         :type other: ``ArithmeticTuple``
1501         """         """
1502         if len(self) != len(other):         if len(self) != len(other):
1503             raise ValueError,"tuple length must match."             raise ValueError("tuple length must match.")
1504         for i in range(len(self)):         for i in range(len(self)):
1505             self.__items[i]-=other[i]             self.__items[i]-=other[i]
1506         return self         return self
# Line 1501  class HomogeneousSaddlePointProblem(obje Line 1528  class HomogeneousSaddlePointProblem(obje
1528        *A* may depend weakly on *v* and *p*.        *A* may depend weakly on *v* and *p*.
1529        """        """
1530        def __init__(self, **kwargs):        def __init__(self, **kwargs):
1531      """          """
1532      initializes the saddle point problem          initializes the saddle point problem
1533      """          """
1534          self.resetControlParameters()          self.resetControlParameters()
1535          self.setTolerance()          self.setTolerance()
1536          self.setAbsoluteTolerance()          self.setAbsoluteTolerance()
1537        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):
1538           """           """
1539           sets a control parameter           sets a control parameter
1540    
1541           :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
1542           :type gamma: ``float``           :type K_p: ``float``
1543           :param gamma_min: minimum value for ``gamma``.           :param K_v: initial value for constant to adjust velocity tolerance
1544           :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``  
1545           :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.
1546           :type rtol_max: ``float``           :type rtol_max: ``float``
1547           :param chi: initial convergence measure.           :param chi_max: maximum tolerable converegence rate.
1548           :type chi: ``float``           :type chi_max: ``float``
1549           :param C_p: initial value for constant to adjust pressure tolerance           :param reduction_factor: reduction factor for adjustment factors.
1550           :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``  
1551           """           """
1552           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)
1553    
1554        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):
1555           """           """
1556           sets a control parameter           sets a control parameter
1557    
1558           :param gamma: ``1/(1-gamma)`` controls the perturbation of the converegence rate due to termination errors in the subproblems.  
1559           :type gamma: ``float``           :param K_p: initial value for constant to adjust pressure tolerance
1560           :param gamma_min: minimum value for ``gamma``.           :type K_p: ``float``
1561           :type gamma_min: ``float``           :param K_v: initial value for constant to adjust velocity tolerance
1562           :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``  
1563           :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.
1564           :type rtol_max: ``float``           :type rtol_max: ``float``
1565           :param chi: initial convergence measure.           :param chi_max: maximum tolerable converegence rate.
1566           :type chi: ``float``           :type chi_max: ``float``
1567           :param C_p: initial value for constant to adjust pressure tolerance           :type reduction_factor: ``float``
1568           :type C_p: ``float``           """
1569           :param C_v: initial value for constant to adjust velocity tolerance           if not K_p == None:
1570           :type C_v: ``float``              if K_p<1:
1571           :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."  
          else:  
             omega_div=self.__omega_div  
   
          if not omega_conv == None:  
             if omega_conv<1:  
                raise ValueError,"omega_conv needs to be greater or equal to 1."  
1572           else:           else:
1573              omega_conv=self.__omega_conv              K_p=self.__K_p
1574    
1575           if not rtol_min == None:           if not K_v == None:
1576              if rtol_min<=0 or rtol_min>=1:              if K_v<1:
1577                 raise ValueError,"rtol_min needs to be positive and less than 1."                 raise ValueError("K_v need to be greater or equal to 1.")
1578           else:           else:
1579              rtol_min=self.__rtol_min              K_v=self.__K_v
1580    
1581           if not rtol_max == None:           if not rtol_max == None:
1582              if rtol_max<=0 or rtol_max>=1:              if rtol_max<=0 or rtol_max>=1:
1583                 raise ValueError,"rtol_max needs to be positive and less than 1."                 raise ValueError("rtol_max needs to be positive and less than 1.")
1584           else:           else:
1585              rtol_max=self.__rtol_max              rtol_max=self.__rtol_max
1586    
1587           if not chi == None:           if not rtol_min == None:
1588              if chi<=0 or chi>1:              if rtol_min<=0 or rtol_min>=1:
1589                 raise ValueError,"chi needs to be positive and less than 1."                 raise ValueError("rtol_min needs to be positive and less than 1.")
1590           else:           else:
1591              chi=self.__chi              rtol_min=self.__rtol_min
1592    
1593           if not C_p == None:           if not chi_max == None:
1594              if C_p<1:              if chi_max<=0 or chi_max>=1:
1595                 raise ValueError,"C_p need to be greater or equal to 1."                 raise ValueError("chi_max needs to be positive and less than 1.")
1596           else:           else:
1597              C_p=self.__C_p              chi_max = self.__chi_max
1598    
1599           if not C_v == None:           if not reduction_factor == None:
1600              if C_v<1:              if reduction_factor<=0 or reduction_factor>1:
1601                 raise ValueError,"C_v need to be greater or equal to 1."                 raise ValueError("reduction_factor need to be between zero and one.")
1602           else:           else:
1603              C_v=self.__C_v              reduction_factor=self.__reduction_factor
1604    
1605           if not safety_factor == None:           if not theta == None:
1606              if safety_factor<=0 or safety_factor>1:              if theta<=0 or theta>1:
1607                 raise ValueError,"safety_factor need to be between zero and one."                 raise ValueError("theta need to be between zero and one.")
1608           else:           else:
1609              safety_factor=self.__safety_factor              theta=self.__theta
1610    
1611           if gamma<gamma_min:           if rtol_min>=rtol_max:
1612                 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  
1613           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  
1614           self.__rtol_max = rtol_max           self.__rtol_max = rtol_max
1615           self.__chi = chi           self.__K_p = K_p
1616           self.__C_p = C_p           self.__K_v = K_v
1617           self.__C_v = C_v           self.__reduction_factor = reduction_factor
1618           self.__safety_factor = safety_factor           self.__theta = theta
1619             self.__rtol_min=rtol_min
1620    
1621        #=============================================================        #=============================================================
1622        def inner_pBv(self,p,Bv):        def inner_pBv(self,p,Bv):
# Line 1657  class HomogeneousSaddlePointProblem(obje Line 1629  class HomogeneousSaddlePointProblem(obje
1629           :rtype: ``float``           :rtype: ``float``
1630           :note: used if PCG is applied.           :note: used if PCG is applied.
1631           """           """
1632           raise NotImplementedError,"no inner product for p and Bv implemented."           raise NotImplementedError("no inner product for p and Bv implemented.")
1633    
1634        def inner_p(self,p0,p1):        def inner_p(self,p0,p1):
1635           """           """
# Line 1668  class HomogeneousSaddlePointProblem(obje Line 1640  class HomogeneousSaddlePointProblem(obje
1640           :return: inner product of p0 and p1           :return: inner product of p0 and p1
1641           :rtype: ``float``           :rtype: ``float``
1642           """           """
1643           raise NotImplementedError,"no inner product for p implemented."           raise NotImplementedError("no inner product for p implemented.")
1644        
1645        def norm_v(self,v):        def norm_v(self,v):
1646           """           """
# Line 1678  class HomogeneousSaddlePointProblem(obje Line 1650  class HomogeneousSaddlePointProblem(obje
1650           :return: norm of v           :return: norm of v
1651           :rtype: non-negative ``float``           :rtype: non-negative ``float``
1652           """           """
1653           raise NotImplementedError,"no norm of v implemented."           raise NotImplementedError("no norm of v implemented.")
1654        def getDV(self, p, v, tol):        def getDV(self, p, v, tol):
1655           """           """
1656           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 1660  class HomogeneousSaddlePointProblem(obje
1660           :return: dv given as *dv= A^{-1} (f-A v-B^*p)*           :return: dv given as *dv= A^{-1} (f-A v-B^*p)*
1661           :note: Only *A* may depend on *v* and *p*           :note: Only *A* may depend on *v* and *p*
1662           """           """
1663           raise NotImplementedError,"no dv calculation implemented."           raise NotImplementedError("no dv calculation implemented.")
1664    
1665                    
1666        def Bv(self,v, tol):        def Bv(self,v, tol):
# Line 1698  class HomogeneousSaddlePointProblem(obje Line 1670  class HomogeneousSaddlePointProblem(obje
1670          :rtype: equal to the type of p          :rtype: equal to the type of p
1671          :note: boundary conditions on p should be zero!          :note: boundary conditions on p should be zero!
1672          """          """
1673          raise NotImplementedError, "no operator B implemented."          raise NotImplementedError("no operator B implemented.")
1674    
1675        def norm_Bv(self,Bv):        def norm_Bv(self,Bv):
1676          """          """
# Line 1707  class HomogeneousSaddlePointProblem(obje Line 1679  class HomogeneousSaddlePointProblem(obje
1679          :rtype: equal to the type of p          :rtype: equal to the type of p
1680          :note: boundary conditions on p should be zero!          :note: boundary conditions on p should be zero!
1681          """          """
1682          raise NotImplementedError, "no norm of Bv implemented."          raise NotImplementedError("no norm of Bv implemented.")
1683    
1684        def solve_AinvBt(self,dp, tol):        def solve_AinvBt(self,dp, tol):
1685           """           """
# Line 1717  class HomogeneousSaddlePointProblem(obje Line 1689  class HomogeneousSaddlePointProblem(obje
1689           :return: the solution of *A dv=B^*dp*           :return: the solution of *A dv=B^*dp*
1690           :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.
1691           """           """
1692           raise NotImplementedError,"no operator A implemented."           raise NotImplementedError("no operator A implemented.")
1693    
1694        def solve_prec(self,Bv, tol):        def solve_prec(self,Bv, tol):
1695           """           """
# Line 1726  class HomogeneousSaddlePointProblem(obje Line 1698  class HomogeneousSaddlePointProblem(obje
1698           :rtype: equal to the type of p           :rtype: equal to the type of p
1699           :note: boundary conditions on p should be zero!           :note: boundary conditions on p should be zero!
1700           """           """
1701           raise NotImplementedError,"no preconditioner for Schur complement implemented."           raise NotImplementedError("no preconditioner for Schur complement implemented.")
1702        #=============================================================        #=============================================================
1703        def __Aprod_PCG(self,dp):        def __Aprod_PCG(self,dp):
1704            dv=self.solve_AinvBt(dp, self.__subtol)            dv=self.solve_AinvBt(dp, self.__subtol)
# Line 1753  class HomogeneousSaddlePointProblem(obje Line 1725  class HomogeneousSaddlePointProblem(obje
1725            :rtype: ``float``            :rtype: ``float``
1726            """            """
1727            f=self.inner_p(p,p)            f=self.inner_p(p,p)
1728            if f<0: raise ValueError,"negative pressure norm."            if f<0: raise ValueError("negative pressure norm.")
1729            return math.sqrt(f)            return math.sqrt(f)
1730                
1731        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 1759  class HomogeneousSaddlePointProblem(obje
1759           self.verbose=verbose           self.verbose=verbose
1760           rtol=self.getTolerance()           rtol=self.getTolerance()
1761           atol=self.getAbsoluteTolerance()           atol=self.getAbsoluteTolerance()
1762    
1763             K_p=self.__K_p
1764             K_v=self.__K_v
1765           correction_step=0           correction_step=0
1766           converged=False           converged=False
          error=None  
1767           chi=None           chi=None
1768           gamma=self.__gamma           eps=None
1769           C_p=self.__C_p  
1770           C_v=self.__C_v           if self.verbose: print(("HomogeneousSaddlePointProblem: start iteration: rtol= %e, atol=%e"%(rtol, atol)))
1771           while not converged:           while not converged:
1772                if error== None or chi == None:  
1773                    gamma_new=gamma/self.__omega_conv               # get tolerance for velecity increment:
1774                else:               if chi == None:
1775                   if chi < self.__chi_max:                  rtol_v=self.__rtol_max
1776                      gamma_new=min(max(gamma*self.__omega_conv,1-chi*error/(self.__safety_factor*ATOL)), 1-chi/self.__chi_max)               else:
1777                   else:                  rtol_v=min(chi/K_v,self.__rtol_max)
1778                      gamma_new=gamma*self.__omega_div               rtol_v=max(rtol_v, self.__rtol_min)
1779                gamma=max(gamma_new, self.__gamma_min)               if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: rtol_v= %e"%(correction_step,rtol_v)))
1780                # calculate velocity for current pressure:               # get velocity increment:
1781                rtol_v=min(max(gamma/(1.+gamma)/C_v,self.__rtol_min), self.__rtol_max)               dv1=self.getDV(p,v,rtol_v)
1782                rtol_p=min(max(gamma/C_p, self.__rtol_min), self.__rtol_max)               v1=v+dv1
1783                self.__subtol=rtol_p**2               Bv1=self.Bv(v1, rtol_v)
1784                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)
1785                if self.verbose: print "HomogeneousSaddlePointProblem: subtolerance: %e"%self.__subtol               norm_dv1=self.norm_v(dv1)
1786                # 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)))
1787                dv1=self.getDV(p,v,rtol_v)               if norm_dv1*self.__theta < norm_Bv1:
1788                v1=v+dv1                  # get tolerance for pressure increment:
1789                Bv1=self.Bv(v1, self.__subtol)                  large_Bv1=True
1790                norm_Bv1=self.norm_Bv(Bv1)                  if chi == None or eps == None:
1791                norm_dv1=self.norm_v(dv1)                     rtol_p=self.__rtol_max
1792                norm_v1=self.norm_v(v1)                  else:
1793                ATOL=norm_v1*rtol+atol                     rtol_p=min(chi**2*eps/K_p/norm_Bv1, self.__rtol_max)
1794                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)
1795                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)))
1796                if max(norm_Bv1, norm_dv1) <= ATOL:                  # now we solve for the pressure increment dp from B*A^{-1}B^* dp = Bv1
1797                    converged=True                  if usePCG:
                   v=v1  
               else:  
                   # now we solve for the pressure increment dp from B*A^{-1}B^* dp = Bv1  
                   if usePCG:  
1798                      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)
1799                      v2=r[0]                      v2=r[0]
1800                      Bv2=r[1]                      Bv2=r[1]
1801                    else:                  else:
1802                        # don't use!!!!
1803                      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)
1804                      dv2=self.solve_AinvBt(dp, self.__subtol)                      dv2=self.solve_AinvBt(dp, self.__subtol)
1805                      v2=v1-dv2                      v2=v1-dv2
1806                      Bv2=self.Bv(v2, self.__subtol)                      Bv2=self.Bv(v2, self.__subtol)
1807                    #                  p2=p+dp
1808                    # convergence indicators:               else:
1809                    #                  large_Bv1=False
1810                    norm_v2=self.norm_v(v2)                  v2=v1
1811                    norm_dv2=self.norm_v(v2-v)                  p2=p
1812                    error_new=max(norm_dv2, norm_Bv1)               # update business:
1813                    norm_Bv2=self.norm_Bv(Bv2)               norm_dv2=self.norm_v(v2-v)
1814                    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)
1815                    if error !=None:               if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: v2 = %e, norm_dv2 = %e"%(correction_step, norm_v2, self.norm_v(v2-v))))
1816                        chi_new=error_new/error               eps, eps_old = max(norm_Bv1, norm_dv2), eps
1817                        if self.verbose: print "HomogeneousSaddlePointProblem: step %s: convergence rate = %e, est. error = %e"%(correction_step,chi_new, error_new)               if eps_old == None:
1818                        if chi != None:                    chi, chi_old = None, chi
1819                            gamma0=max(gamma, 1-chi/chi_new)               else:
1820                            C_p*=gamma0/gamma                    chi, chi_old = min(eps/ eps_old, self.__chi_max), chi
1821                            C_v*=gamma0/gamma*(1+gamma)/(1+gamma0)               if eps != None:
1822                        chi=chi_new                   if chi !=None:
1823                    else:                      if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: convergence rate = %e, correction = %e"%(correction_step,chi, eps)))
1824                        if self.verbose: print "HomogeneousSaddlePointProblem: step %s: est. error = %e"%(correction_step, error_new)                   else:
1825                        if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: correction = %e"%(correction_step, eps)))
1826                    error = error_new               if eps <= rtol*norm_v2+atol :
1827                    correction_step+=1                   converged = True
1828                    if correction_step>max_correction_steps:               else:
1829                          raise CorrectionFailed,"Given up after %d correction steps."%correction_step                   if correction_step>=max_correction_steps:
1830                    v,p=v2,p+dp                        raise CorrectionFailed("Given up after %d correction steps."%correction_step)
1831           if self.verbose: print "HomogeneousSaddlePointProblem: tolerance reached after %s steps."%correction_step                   if chi_old!=None:
1832       return v,p                      K_p=max(1,self.__reduction_factor*K_p,(chi-chi_old)/chi_old**2*K_p)
1833                        K_v=max(1,self.__reduction_factor*K_v,(chi-chi_old)/chi_old**2*K_p)
1834                        if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: new adjustment factor K = %e"%(correction_step,K_p)))
1835                 correction_step+=1
1836                 v,p =v2, p2
1837             if self.verbose: print(("HomogeneousSaddlePointProblem: tolerance reached after %s steps."%correction_step))
1838             return v,p
1839        #========================================================================        #========================================================================
1840        def setTolerance(self,tolerance=1.e-4):        def setTolerance(self,tolerance=1.e-4):
1841           """           """
# Line 1869  class HomogeneousSaddlePointProblem(obje Line 1845  class HomogeneousSaddlePointProblem(obje
1845           :type tolerance: non-negative ``float``           :type tolerance: non-negative ``float``
1846           """           """
1847           if tolerance<0:           if tolerance<0:
1848               raise ValueError,"tolerance must be positive."               raise ValueError("tolerance must be positive.")
1849           self.__rtol=tolerance           self.__rtol=tolerance
1850    
1851        def getTolerance(self):        def getTolerance(self):
# Line 1889  class HomogeneousSaddlePointProblem(obje Line 1865  class HomogeneousSaddlePointProblem(obje
1865           :type tolerance: non-negative ``float``           :type tolerance: non-negative ``float``
1866           """           """
1867           if tolerance<0:           if tolerance<0:
1868               raise ValueError,"tolerance must be non-negative."               raise ValueError("tolerance must be non-negative.")
1869           self.__atol=tolerance           self.__atol=tolerance
1870    
1871        def getAbsoluteTolerance(self):        def getAbsoluteTolerance(self):

Legend:
Removed from v.2793  
changed lines
  Added in v.3981

  ViewVC Help
Powered by ViewVC 1.1.26