/[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 4154 by jfenwick, Tue Jan 22 09:30:23 2013 UTC
# Line 1  Line 1 
1    
2  ########################################################  ##############################################################################
3  #  #
4  # Copyright (c) 2003-2009 by University of Queensland  # Copyright (c) 2003-2013 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-2013 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 escriptcpp
45    escore=escriptcpp
46    from . import linearPDEs
47  import numpy  import numpy
48  import util  from . import util
49  import math  import math
50    
51  class TimeIntegrationManager:  class TimeIntegrationManager(object):
52    """    """
53    A simple mechanism to manage time dependend values.    A simple mechanism to manage time dependend values.
54    
# Line 66  class TimeIntegrationManager: Line 69  class TimeIntegrationManager:
69       Sets up the value manager where ``inital_values`` are the initial values       Sets up the value manager where ``inital_values`` are the initial values
70       and p is the order used for extrapolation.       and p is the order used for extrapolation.
71       """       """
72       if kwargs.has_key("p"):       if "p" in kwargs:
73              self.__p=kwargs["p"]              self.__p=kwargs["p"]
74       else:       else:
75              self.__p=1              self.__p=1
76       if kwargs.has_key("time"):       if "time" in kwargs:
77              self.__t=kwargs["time"]              self.__t=kwargs["time"]
78       else:       else:
79              self.__t=0.              self.__t=0.
# Line 125  class TimeIntegrationManager: Line 128  class TimeIntegrationManager:
128           return out           return out
129    
130    
131  class Projector:  class Projector(object):
132    """    """
133    The Projector is a factory which projects a discontinuous function onto a    The Projector is a factory which projects a discontinuous function onto a
134    continuous function on a given domain.    continuous function on a given domain.
# Line 153  class Projector: Line 156  class Projector:
156      """      """
157      return self.__pde.getSolverOptions()      return self.__pde.getSolverOptions()
158    
159      def getValue(self, input_data):
160        """
161        Projects ``input_data`` onto a continuous function.
162    
163        :param input_data: the data to be projected
164        """
165        return self(input_data)
166    
167    def __call__(self, input_data):    def __call__(self, input_data):
168      """      """
169      Projects ``input_data`` onto a continuous function.      Projects ``input_data`` onto a continuous function.
170    
171      :param input_data: the data to be projected      :param input_data: the data to be projected
172      """      """
173      out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())      out=escore.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())
174      self.__pde.setValue(Y = escript.Data(), Y_reduced = escript.Data())      self.__pde.setValue(Y = escore.Data(), Y_reduced = escore.Data())
175      if input_data.getRank()==0:      if input_data.getRank()==0:
176          self.__pde.setValue(Y = input_data)          self.__pde.setValue(Y = input_data)
177          out=self.__pde.getSolution()          out=self.__pde.getSolution()
# Line 188  class Projector: Line 199  class Projector:
199                      out[i0,i1,i2,i3]=self.__pde.getSolution()                      out[i0,i1,i2,i3]=self.__pde.getSolution()
200      return out      return out
201    
202  class NoPDE:  class NoPDE(object):
203       """       """
204       Solves the following problem for u:       Solves the following problem for u:
205    
# Line 239  class NoPDE: Line 250  class NoPDE:
250           self.__q=q           self.__q=q
251           self.__r=r           self.__r=r
252           self.__u=None           self.__u=None
253           self.__function_space=escript.Solution(self.__domain)           self.__function_space=escore.Solution(self.__domain)
254    
255       def setReducedOn(self):       def setReducedOn(self):
256           """           """
257           Sets the `FunctionSpace` of the solution to `ReducedSolution`.           Sets the `FunctionSpace` of the solution to `ReducedSolution`.
258           """           """
259           self.__function_space=escript.ReducedSolution(self.__domain)           self.__function_space=escore.ReducedSolution(self.__domain)
260           self.__u=None           self.__u=None
261    
262       def setReducedOff(self):       def setReducedOff(self):
263           """           """
264           Sets the `FunctionSpace` of the solution to `Solution`.           Sets the `FunctionSpace` of the solution to `Solution`.
265           """           """
266           self.__function_space=escript.Solution(self.__domain)           self.__function_space=escore.Solution(self.__domain)
267           self.__u=None           self.__u=None
268    
269       def setValue(self,D=None,Y=None,q=None,r=None):       def setValue(self,D=None,Y=None,q=None,r=None):
# Line 291  class NoPDE: Line 302  class NoPDE:
302           """           """
303           if self.__u==None:           if self.__u==None:
304              if self.__D==None:              if self.__D==None:
305                 raise ValueError,"coefficient D is undefined"                 raise ValueError("coefficient D is undefined")
306              D=escript.Data(self.__D,self.__function_space)              D=escore.Data(self.__D,self.__function_space)
307              if D.getRank()>1:              if D.getRank()>1:
308                 raise ValueError,"coefficient D must have rank 0 or 1"                 raise ValueError("coefficient D must have rank 0 or 1")
309              if self.__Y==None:              if self.__Y==None:
310                 self.__u=escript.Data(0.,D.getShape(),self.__function_space)                 self.__u=escore.Data(0.,D.getShape(),self.__function_space)
311              else:              else:
312                 self.__u=util.quotient(self.__Y,D)                 self.__u=1./D*self.__Y
313              if not self.__q==None:              if not self.__q==None:
314                  q=util.wherePositive(escript.Data(self.__q,self.__function_space))                  q=util.wherePositive(escore.Data(self.__q,self.__function_space))
315                  self.__u*=(1.-q)                  self.__u*=(1.-q)
316                  if not self.__r==None: self.__u+=q*self.__r                  if not self.__r==None: self.__u+=q*self.__r
317           return self.__u           return self.__u
318    
319  class Locator:  class Locator(object):
320       """       """
321       Locator provides access to the values of data objects at a given spatial       Locator provides access to the values of data objects at a given spatial
322       coordinate x.       coordinate x.
# Line 325  class Locator: Line 336  class Locator:
336         :param x: location(s) of the Locator         :param x: location(s) of the Locator
337         :type x: ``numpy.ndarray`` or ``list`` of ``numpy.ndarray``         :type x: ``numpy.ndarray`` or ``list`` of ``numpy.ndarray``
338         """         """
339         if isinstance(where,escript.FunctionSpace):         if isinstance(where,escore.FunctionSpace):
340            self.__function_space=where            self.__function_space=where
341         else:         else:
342            self.__function_space=escript.ContinuousFunction(where)            self.__function_space=escore.ContinuousFunction(where)
343         iterative=False         iterative=False
344         if isinstance(x, list):         if isinstance(x, list):
345             if len(x)==0:             if len(x)==0:
346                raise "ValueError", "At least one point must be given."                raise ValueError("At least one point must be given.")
347             try:             try:
348               iter(x[0])               iter(x[0])
349               iterative=True               iterative=True
350             except TypeError:             except TypeError:
351               iterative=False               iterative=False
352           xxx=self.__function_space.getX()
353         if iterative:         if iterative:
354             self.__id=[]             self.__id=[]
355             for p in x:             for p in x:
356                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())
357         else:         else:
358             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()
359    
360       def __str__(self):       def __str__(self):
361         """         """
# Line 400  class Locator: Line 412  class Locator:
412          Returns the value of ``data`` at the Locator if ``data`` is a `Data`          Returns the value of ``data`` at the Locator if ``data`` is a `Data`
413          object otherwise the object is returned.          object otherwise the object is returned.
414          """          """
415          if isinstance(data,escript.Data):          if isinstance(data,escore.Data):
416             dat=util.interpolate(data,self.getFunctionSpace())             dat=util.interpolate(data,self.getFunctionSpace())
417             id=self.getId()             id=self.getId()
418             r=data.getRank()             r=data.getRank()
# Line 421  class Locator: Line 433  class Locator:
433                  return out                  return out
434          else:          else:
435             return data             return data
436              
437         def setValue(self, data, v):
438          """
439          Sets the value of the ``data`` at the Locator.
440          """
441          if isinstance(data, escore.Data):
442             if data.getFunctionSpace()!=self.getFunctionSpace():
443               raise TypeError("setValue: FunctionSpace of Locator and Data object must match.")
444             data.expand()  
445             id=self.getId()
446             if isinstance(id, list):
447              for i in id:
448               data._setTupleForGlobalDataPoint(i[1], i[0], v)
449             else:
450               data._setTupleForGlobalDataPoint(id[1], id[0], v)
451          else:
452               raise TypeError("setValue: Invalid argument type.")
453    
454    
455  def getInfLocator(arg):  def getInfLocator(arg):
456      """      """
457      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.
458      """      """
459      if not isinstance(arg, escript.Data):      if not isinstance(arg, escore.Data):
460      raise TypeError,"getInfLocator: Unknown argument type."         raise TypeError("getInfLocator: Unknown argument type.")
461      a_inf=util.inf(arg)      a_inf=util.inf(arg)
462      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
463      x=arg.getFunctionSpace().getX()      x=arg.getFunctionSpace().getX()
# Line 439  def getSupLocator(arg): Line 468  def getSupLocator(arg):
468      """      """
469      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.
470      """      """
471      if not isinstance(arg, escript.Data):      if not isinstance(arg, escore.Data):
472      raise TypeError,"getInfLocator: Unknown argument type."         raise TypeError("getInfLocator: Unknown argument type.")
473      a_inf=util.sup(arg)      a_inf=util.sup(arg)
474      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
475      x=arg.getFunctionSpace().getX()      x=arg.getFunctionSpace().getX()
# Line 539  def PCG(r, Aprod, x, Msolve, bilinearfor Line 568  def PCG(r, Aprod, x, Msolve, bilinearfor
568     rhat=Msolve(r)     rhat=Msolve(r)
569     d = rhat     d = rhat
570     rhat_dot_r = bilinearform(rhat, r)     rhat_dot_r = bilinearform(rhat, r)
571     if rhat_dot_r<0: raise NegativeNorm,"negative norm."     if rhat_dot_r<0: raise NegativeNorm("negative norm.")
572     norm_r0=math.sqrt(rhat_dot_r)     norm_r0=math.sqrt(rhat_dot_r)
573     atol2=atol+rtol*norm_r0     atol2=atol+rtol*norm_r0
574     if atol2<=0:     if atol2<=0:
575        raise ValueError,"Non-positive tolarance."        raise ValueError("Non-positive tolarance.")
576     atol2=max(atol2, 100. * util.EPSILON * norm_r0)     atol2=max(atol2, 100. * util.EPSILON * norm_r0)
577    
578     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)))
579    
580    
581     while not math.sqrt(rhat_dot_r) <= atol2:     while not math.sqrt(rhat_dot_r) <= atol2:
582         iter+=1         iter+=1
583         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)
584    
585         q=Aprod(d)         q=Aprod(d)
586         alpha = rhat_dot_r / bilinearform(d, q)         alpha = rhat_dot_r / bilinearform(d, q)
587         x += alpha * d         x += alpha * d
588         if isinstance(q,ArithmeticTuple):         if isinstance(q,ArithmeticTuple):
589         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__
590         else:         else:
591             r += (-alpha) * q             r += (-alpha) * q
592         rhat=Msolve(r)         rhat=Msolve(r)
# Line 567  def PCG(r, Aprod, x, Msolve, bilinearfor Line 596  def PCG(r, Aprod, x, Msolve, bilinearfor
596         d=rhat         d=rhat
597    
598         rhat_dot_r = rhat_dot_r_new         rhat_dot_r = rhat_dot_r_new
599         if rhat_dot_r<0: raise NegativeNorm,"negative norm."         if rhat_dot_r<0: raise NegativeNorm("negative norm.")
600         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))))
601     if verbose: print "PCG: tolerance reached after %s steps."%iter     if verbose: print(("PCG: tolerance reached after %s steps."%iter))
602     return x,r,math.sqrt(rhat_dot_r)     return x,r,math.sqrt(rhat_dot_r)
603    
604  class Defect(object):  class Defect(object):
# Line 603  class Defect(object): Line 632  class Defect(object):
632          :note: by default ``sqrt(self.bilinearform(x,x)`` is returned.          :note: by default ``sqrt(self.bilinearform(x,x)`` is returned.
633          """          """
634          s=self.bilinearform(x,x)          s=self.bilinearform(x,x)
635          if s<0: raise NegativeNorm,"negative norm."          if s<0: raise NegativeNorm("negative norm.")
636          return math.sqrt(s)          return math.sqrt(s)
637    
638      def eval(self,x):      def eval(self,x):
# Line 627  class Defect(object): Line 656  class Defect(object):
656          :param inc: relative increment length          :param inc: relative increment length
657          :type inc: positive ``float``          :type inc: positive ``float``
658          """          """
659          if inc<=0: raise ValueError,"positive increment required."          if inc<=0: raise ValueError("positive increment required.")
660          self.__inc=inc          self.__inc=inc
661    
662      def getDerivativeIncrementLength(self):      def getDerivativeIncrementLength(self):
# Line 698  def NewtonGMRES(defect, x, iter_max=100, Line 727  def NewtonGMRES(defect, x, iter_max=100,
727     :rtype: same type as the initial guess     :rtype: same type as the initial guess
728     """     """
729     lmaxit=iter_max     lmaxit=iter_max
730     if atol<0: raise ValueError,"atol needs to be non-negative."     if atol<0: raise ValueError("atol needs to be non-negative.")
731     if rtol<0: raise ValueError,"rtol needs to be non-negative."     if rtol<0: raise ValueError("rtol needs to be non-negative.")
732     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.")
733     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)
734     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)
735    
736     F=defect(x)     F=defect(x)
737     fnrm=defect.norm(F)     fnrm=defect.norm(F)
738     stop_tol=atol + rtol*fnrm     stop_tol=atol + rtol*fnrm
739     subtol=subtol_max     subtol=subtol_max
740     if verbose: print "NewtonGMRES: initial residual = %e."%fnrm     if verbose: print(("NewtonGMRES: initial residual = %e."%fnrm))
741     if verbose: print "             tolerance = %e."%subtol     if verbose: print(("             tolerance = %e."%subtol))
742     iter=1     iter=1
743     #     #
744     # main iteration loop     # main iteration loop
745     #     #
746     while not fnrm<=stop_tol:     while not fnrm<=stop_tol:
747              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)
748              #              #
749          #   adjust subtol_          #   adjust subtol_
750          #          #
751              if iter > 1:              if iter > 1:
752             rat=fnrm/fnrmo                 rat=fnrm/fnrmo
753                 subtol_old=subtol                 subtol_old=subtol
754             subtol=gamma*rat**2                 subtol=gamma*rat**2
755             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)
756             subtol=max(min(subtol,subtol_max), .5*stop_tol/fnrm)                 subtol=max(min(subtol,subtol_max), .5*stop_tol/fnrm)
757          #          #
758          # calculate newton increment xc          # calculate newton increment xc
759              #     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 761  def NewtonGMRES(defect, x, iter_max=100,
761              #     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
762              #              #
763              #              #
764              if verbose: print "             subiteration (GMRES) is called with relative tolerance %e."%subtol              if verbose: print(("             subiteration (GMRES) is called with relative tolerance %e."%subtol))
765              try:              try:
766                 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)
767              except MaxIterReached:              except MaxIterReached:
768                 raise MaxIterReached,"maximum number of %s steps reached."%iter_max                 raise MaxIterReached("maximum number of %s steps reached."%iter_max)
769              if sub_iter<0:              if sub_iter<0:
770                 iter+=sub_iter_max                 iter+=sub_iter_max
771              else:              else:
772                 iter+=sub_iter                 iter+=sub_iter
773              # ====              # ====
774          x+=xc              x+=xc
775              F=defect(x)              F=defect(x)
776          iter+=1              iter+=1
777              fnrmo, fnrm=fnrm, defect.norm(F)              fnrmo, fnrm=fnrm, defect.norm(F)
778              if verbose: print "             step %s: residual %e."%(iter,fnrm)              if verbose: print(("             step %s: residual %e."%(iter,fnrm)))
779     if verbose: print "NewtonGMRES: completed after %s steps."%iter     if verbose: print(("NewtonGMRES: completed after %s steps."%iter))
780     return x     return x
781    
782  def __givapp(c,s,vin):  def __givapp(c,s,vin):
# Line 763  def __givapp(c,s,vin): Line 792  def __givapp(c,s,vin):
792      else:      else:
793          for i in range(len(c)):          for i in range(len(c)):
794              w1=c[i]*vrot[i]-s[i]*vrot[i+1]              w1=c[i]*vrot[i]-s[i]*vrot[i+1]
795          w2=s[i]*vrot[i]+c[i]*vrot[i+1]              w2=s[i]*vrot[i]+c[i]*vrot[i+1]
796              vrot[i]=w1              vrot[i]=w1
797              vrot[i+1]=w2              vrot[i+1]=w2
798      return vrot      return vrot
# Line 783  def __FDGMRES(F0, defect, x0, atol, iter Line 812  def __FDGMRES(F0, defect, x0, atol, iter
812     iter=0     iter=0
813     while rho > atol and iter<iter_restart-1:     while rho > atol and iter<iter_restart-1:
814          if iter  >= iter_max:          if iter  >= iter_max:
815              raise MaxIterReached,"maximum number of %s steps reached."%iter_max              raise MaxIterReached("maximum number of %s steps reached."%iter_max)
816    
817          p=defect.derivative(F0,x0,v[iter], v_is_normalised=True)          p=defect.derivative(F0,x0,v[iter], v_is_normalised=True)
818          v.append(p)          v.append(p)
# Line 845  def __FDGMRES(F0, defect, x0, atol, iter Line 874  def __FDGMRES(F0, defect, x0, atol, iter
874            i=i-1            i=i-1
875       xhat=v[iter-1]*y[iter-1]       xhat=v[iter-1]*y[iter-1]
876       for i in range(iter-1):       for i in range(iter-1):
877      xhat += v[i]*y[i]         xhat += v[i]*y[i]
878     else :     else :
879        xhat=v[0] * 0        xhat=v[0] * 0
880    
# Line 904  def GMRES(r, Aprod, x, bilinearform, ato Line 933  def GMRES(r, Aprod, x, bilinearform, ato
933     iter=0     iter=0
934     if rtol>0:     if rtol>0:
935        r_dot_r = bilinearform(r, r)        r_dot_r = bilinearform(r, r)
936        if r_dot_r<0: raise NegativeNorm,"negative norm."        if r_dot_r<0: raise NegativeNorm("negative norm.")
937        atol2=atol+rtol*math.sqrt(r_dot_r)        atol2=atol+rtol*math.sqrt(r_dot_r)
938        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)))
939     else:     else:
940        atol2=atol        atol2=atol
941        if verbose: print "GMRES: absolute tolerance = %e"%atol2        if verbose: print(("GMRES: absolute tolerance = %e"%atol2))
942     if atol2<=0:     if atol2<=0:
943        raise ValueError,"Non-positive tolarance."        raise ValueError("Non-positive tolarance.")
944    
945     while True:     while True:
946        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)
947        if restarted:        if restarted:
948           r2 = r-Aprod(x-x2)           r2 = r-Aprod(x-x2)
949        else:        else:
# Line 923  def GMRES(r, Aprod, x, bilinearform, ato Line 952  def GMRES(r, Aprod, x, bilinearform, ato
952        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)
953        iter+=iter_restart        iter+=iter_restart
954        if stopped: break        if stopped: break
955        if verbose: print "GMRES: restart."        if verbose: print("GMRES: restart.")
956        restarted=True        restarted=True
957     if verbose: print "GMRES: tolerance has been reached."     if verbose: print("GMRES: tolerance has been reached.")
958     return x     return x
959    
960  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 967  def _GMRESm(r, Aprod, x, bilinearform, a
967     v=[]     v=[]
968    
969     r_dot_r = bilinearform(r, r)     r_dot_r = bilinearform(r, r)
970     if r_dot_r<0: raise NegativeNorm,"negative norm."     if r_dot_r<0: raise NegativeNorm("negative norm.")
971     rho=math.sqrt(r_dot_r)     rho=math.sqrt(r_dot_r)
972    
973     v.append(r/rho)     v.append(r/rho)
974     g[0]=rho     g[0]=rho
975    
976     if verbose: print "GMRES: initial residual %e (absolute tolerance = %e)"%(rho,atol)     if verbose: print(("GMRES: initial residual %e (absolute tolerance = %e)"%(rho,atol)))
977     while not (rho<=atol or iter==iter_restart):     while not (rho<=atol or iter==iter_restart):
978    
979      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)
980    
981          if P_R!=None:          if P_R!=None:
982              p=Aprod(P_R(v[iter]))              p=Aprod(P_R(v[iter]))
983          else:          else:
984          p=Aprod(v[iter])              p=Aprod(v[iter])
985      v.append(p)          v.append(p)
986    
987      v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))          v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
988    
989  # Modified Gram-Schmidt  # Modified Gram-Schmidt
990      for j in range(iter+1):          for j in range(iter+1):
991        h[j,iter]=bilinearform(v[j],v[iter+1])            h[j,iter]=bilinearform(v[j],v[iter+1])
992        v[iter+1]-=h[j,iter]*v[j]            v[iter+1]-=h[j,iter]*v[j]
993    
994      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]))
995      v_norm2=h[iter+1,iter]          v_norm2=h[iter+1,iter]
996    
997  # Reorthogonalize if needed  # Reorthogonalize if needed
998      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)
999       for j in range(iter+1):           for j in range(iter+1):
1000          hr=bilinearform(v[j],v[iter+1])              hr=bilinearform(v[j],v[iter+1])
1001              h[j,iter]=h[j,iter]+hr              h[j,iter]=h[j,iter]+hr
1002              v[iter+1] -= hr*v[j]              v[iter+1] -= hr*v[j]
1003    
1004       v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))           v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
1005       h[iter+1,iter]=v_norm2           h[iter+1,iter]=v_norm2
1006    
1007  #   watch out for happy breakdown  #   watch out for happy breakdown
1008          if not v_norm2 == 0:          if not v_norm2 == 0:
1009           v[iter+1]=v[iter+1]/h[iter+1,iter]           v[iter+1]=v[iter+1]/h[iter+1,iter]
1010    
1011  #   Form and store the information for the new Givens rotation  #   Form and store the information for the new Givens rotation
1012      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])
1013      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])
1014    
1015      if mu!=0 :          if mu!=0 :
1016          c[iter]=h[iter,iter]/mu                  c[iter]=h[iter,iter]/mu
1017          s[iter]=-h[iter+1,iter]/mu                  s[iter]=-h[iter+1,iter]/mu
1018          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]
1019          h[iter+1,iter]=0.0                  h[iter+1,iter]=0.0
1020                  gg=__givapp(c[iter],s[iter],[g[iter],g[iter+1]])                  gg=__givapp(c[iter],s[iter],[g[iter],g[iter+1]])
1021                  g[iter]=gg[0]                  g[iter]=gg[0]
1022                  g[iter+1]=gg[1]                  g[iter+1]=gg[1]
1023  # Update the residual norm  # Update the residual norm
1024    
1025          rho=abs(g[iter+1])          rho=abs(g[iter+1])
1026          if verbose: print "GMRES: iteration step %s: residual %e"%(iter,rho)          if verbose: print(("GMRES: iteration step %s: residual %e"%(iter,rho)))
1027      iter+=1          iter+=1
1028    
1029  # At this point either iter > iter_max or rho < tol.  # At this point either iter > iter_max or rho < tol.
1030  # It's time to compute x and leave.  # It's time to compute x and leave.
1031    
1032     if verbose: print "GMRES: iteration stopped after %s step."%iter     if verbose: print(("GMRES: iteration stopped after %s step."%iter))
1033     if iter > 0 :     if iter > 0 :
1034       y=numpy.zeros(iter,numpy.float64)       y=numpy.zeros(iter,numpy.float64)
1035       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 1040  def _GMRESm(r, Aprod, x, bilinearform, a
1040            i=i-1            i=i-1
1041       xhat=v[iter-1]*y[iter-1]       xhat=v[iter-1]*y[iter-1]
1042       for i in range(iter-1):       for i in range(iter-1):
1043      xhat += v[i]*y[i]         xhat += v[i]*y[i]
1044     else:     else:
1045       xhat=v[0] * 0       xhat=v[0] * 0
1046     if P_R!=None:     if P_R!=None:
# Line 1084  def MINRES(r, Aprod, x, Msolve, bilinear Line 1113  def MINRES(r, Aprod, x, Msolve, bilinear
1113      y = Msolve(r)      y = Msolve(r)
1114      beta1 = bilinearform(y,r)      beta1 = bilinearform(y,r)
1115    
1116      if beta1< 0: raise NegativeNorm,"negative norm."      if beta1< 0: raise NegativeNorm("negative norm.")
1117    
1118      #  If r = 0 exactly, stop with x      #  If r = 0 exactly, stop with x
1119      if beta1==0: return x      if beta1==0: return x
# Line 1119  def MINRES(r, Aprod, x, Msolve, bilinear Line 1148  def MINRES(r, Aprod, x, Msolve, bilinear
1148      # --------------------------------------------------------------------      # --------------------------------------------------------------------
1149      while not rnorm<=atol+rtol*Anorm*ynorm:    #  checks ||r|| < (||A|| ||x||) * TOL      while not rnorm<=atol+rtol*Anorm*ynorm:    #  checks ||r|| < (||A|| ||x||) * TOL
1150    
1151      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)
1152          iter    = iter  +  1          iter    = iter  +  1
1153    
1154          #-----------------------------------------------------------------          #-----------------------------------------------------------------
# Line 1149  def MINRES(r, Aprod, x, Msolve, bilinear Line 1178  def MINRES(r, Aprod, x, Msolve, bilinear
1178          y = Msolve(r2)          y = Msolve(r2)
1179          oldb   = beta                           # oldb = betak          oldb   = beta                           # oldb = betak
1180          beta   = bilinearform(y,r2)             # beta = betak+1^2          beta   = bilinearform(y,r2)             # beta = betak+1^2
1181          if beta < 0: raise NegativeNorm,"negative norm."          if beta < 0: raise NegativeNorm("negative norm.")
1182    
1183          beta   = math.sqrt( beta )          beta   = math.sqrt( beta )
1184          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 1285  def TFQMR(r, Aprod, x, bilinearform, ato
1285    theta = 0.0;    theta = 0.0;
1286    eta = 0.0;    eta = 0.0;
1287    rho=bilinearform(r,r)    rho=bilinearform(r,r)
1288    if rho < 0: raise NegativeNorm,"negative norm."    if rho < 0: raise NegativeNorm("negative norm.")
1289    tau = math.sqrt(rho)    tau = math.sqrt(rho)
1290    norm_r0=tau    norm_r0=tau
1291    while tau>atol+rtol*norm_r0:    while tau>atol+rtol*norm_r0:
1292      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)
1293    
1294      sigma = bilinearform(r,v)      sigma = bilinearform(r,v)
1295      if sigma == 0.0: raise IterationBreakDown,'TFQMR breakdown, sigma=0'      if sigma == 0.0: raise IterationBreakDown('TFQMR breakdown, sigma=0')
1296    
1297      alpha = rho / sigma      alpha = rho / sigma
1298    
# Line 1291  def TFQMR(r, Aprod, x, bilinearform, ato Line 1320  def TFQMR(r, Aprod, x, bilinearform, ato
1320  #  #
1321  #  Try to terminate the iteration at each pass through the loop  #  Try to terminate the iteration at each pass through the loop
1322  #  #
1323      if rho == 0.0: raise IterationBreakDown,'TFQMR breakdown, rho=0'      if rho == 0.0: raise IterationBreakDown('TFQMR breakdown, rho=0')
1324    
1325      rhon = bilinearform(r,w)      rhon = bilinearform(r,w)
1326      beta = rhon / rho;      beta = rhon / rho;
# Line 1316  class ArithmeticTuple(object): Line 1345  class ArithmeticTuple(object):
1345    
1346         from esys.escript import Data         from esys.escript import Data
1347         from numpy import array         from numpy import array
1348         a=Data(...)         a=eData(...)
1349         b=array([1.,4.])         b=array([1.,4.])
1350         x=ArithmeticTuple(a,b)         x=ArithmeticTuple(a,b)
1351         y=5.*x         y=5.*x
# Line 1363  class ArithmeticTuple(object): Line 1392  class ArithmeticTuple(object):
1392         try:         try:
1393             l=len(other)             l=len(other)
1394             if l!=len(self):             if l!=len(self):
1395                 raise ValueError,"length of arguments don't match."                 raise ValueError("length of arguments don't match.")
1396             for i in range(l): out.append(self[i]*other[i])             for i in range(l):
1397            if self.__isEmpty(self[i]) or self.__isEmpty(other[i]):
1398                out.append(escore.Data())
1399            else:
1400                out.append(self[i]*other[i])
1401         except TypeError:         except TypeError:
1402             for i in range(len(self)): out.append(self[i]*other)          for i in range(len(self)):  
1403            if self.__isEmpty(self[i]) or self.__isEmpty(other):
1404                out.append(escore.Data())
1405            else:
1406                out.append(self[i]*other)
1407         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1408    
1409     def __rmul__(self,other):     def __rmul__(self,other):
1410         """        """
1411         Scales by ``other`` from the left.        Scales by ``other`` from the left.
1412    
1413         :param other: scaling factor        :param other: scaling factor
1414         :type other: ``float``        :type other: ``float``
1415         :return: itemwise other*self        :return: itemwise other*self
1416         :rtype: `ArithmeticTuple`        :rtype: `ArithmeticTuple`
1417         """        """
1418         out=[]        out=[]
1419         try:        try:
1420             l=len(other)        l=len(other)
1421             if l!=len(self):        if l!=len(self):
1422                 raise ValueError,"length of arguments don't match."            raise ValueError("length of arguments don't match.")
1423             for i in range(l): out.append(other[i]*self[i])        for i in range(l):
1424         except TypeError:          if self.__isEmpty(self[i]) or self.__isEmpty(other[i]):
1425             for i in range(len(self)): out.append(other*self[i])              out.append(escore.Data())
1426         return ArithmeticTuple(*tuple(out))          else:
1427                out.append(other[i]*self[i])
1428          except TypeError:
1429          for i in range(len(self)):  
1430            if self.__isEmpty(self[i]) or self.__isEmpty(other):
1431                out.append(escore.Data())
1432            else:
1433                out.append(other*self[i])
1434          return ArithmeticTuple(*tuple(out))
1435    
1436     def __div__(self,other):     def __div__(self,other):
1437         """         """
# Line 1400  class ArithmeticTuple(object): Line 1445  class ArithmeticTuple(object):
1445         return self*(1/other)         return self*(1/other)
1446    
1447     def __rdiv__(self,other):     def __rdiv__(self,other):
1448         """        """
1449         Scales by (1/``other``) from the left.        Scales by (1/``other``) from the left.
1450    
1451         :param other: scaling factor        :param other: scaling factor
1452         :type other: ``float``        :type other: ``float``
1453         :return: itemwise other/self        :return: itemwise other/self
1454         :rtype: `ArithmeticTuple`        :rtype: `ArithmeticTuple`
1455         """        """
1456         out=[]        out=[]
1457         try:        try:
1458             l=len(other)        l=len(other)
1459             if l!=len(self):        if l!=len(self):
1460                 raise ValueError,"length of arguments don't match."            raise ValueError("length of arguments don't match.")
1461             for i in range(l): out.append(other[i]/self[i])        
1462         except TypeError:        for i in range(l):
1463             for i in range(len(self)): out.append(other/self[i])          if self.__isEmpty(self[i]):
1464         return ArithmeticTuple(*tuple(out))              raise ZeroDivisionError("in component %s"%i)
1465            else:
1466                if self.__isEmpty(other[i]):
1467                out.append(escore.Data())
1468                else:
1469                out.append(other[i]/self[i])
1470          except TypeError:
1471          for i in range(len(self)):
1472            if self.__isEmpty(self[i]):
1473                raise ZeroDivisionError("in component %s"%i)
1474            else:
1475                if self.__isEmpty(other):
1476                out.append(escore.Data())
1477                else:
1478                out.append(other/self[i])
1479          return ArithmeticTuple(*tuple(out))
1480    
1481     def __iadd__(self,other):     def __iadd__(self,other):
1482         """        """
1483         Inplace addition of ``other`` to self.        Inplace addition of ``other`` to self.
1484    
1485         :param other: increment        :param other: increment
1486         :type other: ``ArithmeticTuple``        :type other: ``ArithmeticTuple``
1487         """        """
1488         if len(self) != len(other):        if len(self) != len(other):
1489             raise ValueError,"tuple lengths must match."        raise ValueError("tuple lengths must match.")
1490         for i in range(len(self)):        for i in range(len(self)):
1491             self.__items[i]+=other[i]        if self.__isEmpty(self.__items[i]):
1492         return self            self.__items[i]=other[i]
1493          else:
1494              self.__items[i]+=other[i]
1495              
1496          return self
1497    
1498     def __add__(self,other):     def __add__(self,other):
1499         """        """
1500         Adds ``other`` to self.        Adds ``other`` to self.
1501    
1502         :param other: increment        :param other: increment
1503         :type other: ``ArithmeticTuple``        :type other: ``ArithmeticTuple``
1504         """        """
1505         out=[]        out=[]
1506         try:        try:
1507             l=len(other)        l=len(other)
1508             if l!=len(self):        if l!=len(self):
1509                 raise ValueError,"length of arguments don't match."            raise ValueError("length of arguments don't match.")
1510             for i in range(l): out.append(self[i]+other[i])        for i in range(l):
1511         except TypeError:          if self.__isEmpty(self[i]):
1512             for i in range(len(self)): out.append(self[i]+other)              out.append(other[i])
1513         return ArithmeticTuple(*tuple(out))          elif self.__isEmpty(other[i]):
1514                out.append(self[i])
1515            else:
1516                out.append(self[i]+other[i])
1517          except TypeError:
1518            for i in range(len(self)):    
1519            if self.__isEmpty(self[i]):
1520                out.append(other)
1521            elif self.__isEmpty(other):
1522                out.append(self[i])
1523            else:
1524                out.append(self[i]+other)
1525          return ArithmeticTuple(*tuple(out))
1526    
1527     def __sub__(self,other):     def __sub__(self,other):
1528         """        """
1529         Subtracts ``other`` from self.        Subtracts ``other`` from self.
1530    
1531         :param other: decrement        :param other: decrement
1532         :type other: ``ArithmeticTuple``        :type other: ``ArithmeticTuple``
1533         """        """
1534         out=[]        out=[]
1535         try:        try:
1536             l=len(other)        l=len(other)
1537             if l!=len(self):        if l!=len(self):
1538                 raise ValueError,"length of arguments don't match."            raise ValueError("length of arguments don't match.")
1539             for i in range(l): out.append(self[i]-other[i])        for i in range(l):
1540         except TypeError:          if self.__isEmpty(other[i]):
1541             for i in range(len(self)): out.append(self[i]-other)              out.append(self[i])
1542         return ArithmeticTuple(*tuple(out))          elif self.__isEmpty(self[i]):
1543                out.append(-other[i])
1544            else:
1545                out.append(self[i]-other[i])
1546          except TypeError:
1547            for i in range(len(self)):    
1548            if  self.__isEmpty(other):
1549                out.append(self[i])
1550            elif self.__isEmpty(self[i]):
1551                out.append(-other)
1552            else:
1553                out.append(self[i]-other)
1554                
1555          return ArithmeticTuple(*tuple(out))
1556    
1557     def __isub__(self,other):     def __isub__(self,other):
1558         """        """
1559         Inplace subtraction of ``other`` from self.        Inplace subtraction of ``other`` from self.
1560    
1561         :param other: decrement        :param other: decrement
1562         :type other: ``ArithmeticTuple``        :type other: ``ArithmeticTuple``
1563         """        """
1564         if len(self) != len(other):        if len(self) != len(other):
1565             raise ValueError,"tuple length must match."        raise ValueError("tuple length must match.")
1566         for i in range(len(self)):        for i in range(len(self)):
1567             self.__items[i]-=other[i]        if not self.__isEmpty(other[i]):
1568         return self            if self.__isEmpty(self.__items[i]):
1569              self.__items[i]=-other[i]
1570              else:
1571              self.__items[i]=other[i]
1572          return self
1573    
1574     def __neg__(self):     def __neg__(self):
1575         """        """
1576         Negates values.        Negates values.
1577         """        """
1578         out=[]        out=[]
1579         for i in range(len(self)):        for i in range(len(self)):
1580             out.append(-self[i])        if self.__isEmpty(self[i]):
1581         return ArithmeticTuple(*tuple(out))            out.append(escore.Data())
1582          else:
1583              out.append(-self[i])
1584          
1585          return ArithmeticTuple(*tuple(out))
1586       def __isEmpty(self, d):
1587        if isinstance(d, escore.Data):
1588        return d.isEmpty()
1589        else:
1590        return False
1591    
1592    
1593  class HomogeneousSaddlePointProblem(object):  class HomogeneousSaddlePointProblem(object):
# Line 1501  class HomogeneousSaddlePointProblem(obje Line 1603  class HomogeneousSaddlePointProblem(obje
1603        *A* may depend weakly on *v* and *p*.        *A* may depend weakly on *v* and *p*.
1604        """        """
1605        def __init__(self, **kwargs):        def __init__(self, **kwargs):
1606      """          """
1607      initializes the saddle point problem          initializes the saddle point problem
1608      """          """
1609          self.resetControlParameters()          self.resetControlParameters()
1610          self.setTolerance()          self.setTolerance()
1611          self.setAbsoluteTolerance()          self.setAbsoluteTolerance()
1612        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):
1613           """           """
1614           sets a control parameter           sets a control parameter
1615    
1616           :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
1617           :type gamma: ``float``           :type K_p: ``float``
1618           :param gamma_min: minimum value for ``gamma``.           :param K_v: initial value for constant to adjust velocity tolerance
1619           :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``  
1620           :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.
1621           :type rtol_max: ``float``           :type rtol_max: ``float``
1622           :param chi: initial convergence measure.           :param chi_max: maximum tolerable converegence rate.
1623           :type chi: ``float``           :type chi_max: ``float``
1624           :param C_p: initial value for constant to adjust pressure tolerance           :param reduction_factor: reduction factor for adjustment factors.
1625           :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``  
1626           """           """
1627           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)
1628    
1629        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):
1630           """           """
1631           sets a control parameter           sets a control parameter
1632    
1633           :param gamma: ``1/(1-gamma)`` controls the perturbation of the converegence rate due to termination errors in the subproblems.  
1634           :type gamma: ``float``           :param K_p: initial value for constant to adjust pressure tolerance
1635           :param gamma_min: minimum value for ``gamma``.           :type K_p: ``float``
1636           :type gamma_min: ``float``           :param K_v: initial value for constant to adjust velocity tolerance
1637           :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``  
1638           :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.
1639           :type rtol_max: ``float``           :type rtol_max: ``float``
1640           :param chi: initial convergence measure.           :param chi_max: maximum tolerable converegence rate.
1641           :type chi: ``float``           :type chi_max: ``float``
1642           :param C_p: initial value for constant to adjust pressure tolerance           :type reduction_factor: ``float``
1643           :type C_p: ``float``           """
1644           :param C_v: initial value for constant to adjust velocity tolerance           if not K_p == None:
1645           :type C_v: ``float``              if K_p<1:
1646           :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."  
1647           else:           else:
1648              omega_conv=self.__omega_conv              K_p=self.__K_p
1649    
1650           if not rtol_min == None:           if not K_v == None:
1651              if rtol_min<=0 or rtol_min>=1:              if K_v<1:
1652                 raise ValueError,"rtol_min needs to be positive and less than 1."                 raise ValueError("K_v need to be greater or equal to 1.")
1653           else:           else:
1654              rtol_min=self.__rtol_min              K_v=self.__K_v
1655    
1656           if not rtol_max == None:           if not rtol_max == None:
1657              if rtol_max<=0 or rtol_max>=1:              if rtol_max<=0 or rtol_max>=1:
1658                 raise ValueError,"rtol_max needs to be positive and less than 1."                 raise ValueError("rtol_max needs to be positive and less than 1.")
1659           else:           else:
1660              rtol_max=self.__rtol_max              rtol_max=self.__rtol_max
1661    
1662           if not chi == None:           if not rtol_min == None:
1663              if chi<=0 or chi>1:              if rtol_min<=0 or rtol_min>=1:
1664                 raise ValueError,"chi needs to be positive and less than 1."                 raise ValueError("rtol_min needs to be positive and less than 1.")
1665           else:           else:
1666              chi=self.__chi              rtol_min=self.__rtol_min
1667    
1668           if not C_p == None:           if not chi_max == None:
1669              if C_p<1:              if chi_max<=0 or chi_max>=1:
1670                 raise ValueError,"C_p need to be greater or equal to 1."                 raise ValueError("chi_max needs to be positive and less than 1.")
1671           else:           else:
1672              C_p=self.__C_p              chi_max = self.__chi_max
1673    
1674           if not C_v == None:           if not reduction_factor == None:
1675              if C_v<1:              if reduction_factor<=0 or reduction_factor>1:
1676                 raise ValueError,"C_v need to be greater or equal to 1."                 raise ValueError("reduction_factor need to be between zero and one.")
1677           else:           else:
1678              C_v=self.__C_v              reduction_factor=self.__reduction_factor
1679    
1680           if not safety_factor == None:           if not theta == None:
1681              if safety_factor<=0 or safety_factor>1:              if theta<=0 or theta>1:
1682                 raise ValueError,"safety_factor need to be between zero and one."                 raise ValueError("theta need to be between zero and one.")
1683           else:           else:
1684              safety_factor=self.__safety_factor              theta=self.__theta
1685    
1686           if gamma<gamma_min:           if rtol_min>=rtol_max:
1687                 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  
1688           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  
1689           self.__rtol_max = rtol_max           self.__rtol_max = rtol_max
1690           self.__chi = chi           self.__K_p = K_p
1691           self.__C_p = C_p           self.__K_v = K_v
1692           self.__C_v = C_v           self.__reduction_factor = reduction_factor
1693           self.__safety_factor = safety_factor           self.__theta = theta
1694             self.__rtol_min=rtol_min
1695    
1696        #=============================================================        #=============================================================
1697        def inner_pBv(self,p,Bv):        def inner_pBv(self,p,Bv):
# Line 1657  class HomogeneousSaddlePointProblem(obje Line 1704  class HomogeneousSaddlePointProblem(obje
1704           :rtype: ``float``           :rtype: ``float``
1705           :note: used if PCG is applied.           :note: used if PCG is applied.
1706           """           """
1707           raise NotImplementedError,"no inner product for p and Bv implemented."           raise NotImplementedError("no inner product for p and Bv implemented.")
1708    
1709        def inner_p(self,p0,p1):        def inner_p(self,p0,p1):
1710           """           """
# Line 1668  class HomogeneousSaddlePointProblem(obje Line 1715  class HomogeneousSaddlePointProblem(obje
1715           :return: inner product of p0 and p1           :return: inner product of p0 and p1
1716           :rtype: ``float``           :rtype: ``float``
1717           """           """
1718           raise NotImplementedError,"no inner product for p implemented."           raise NotImplementedError("no inner product for p implemented.")
1719        
1720        def norm_v(self,v):        def norm_v(self,v):
1721           """           """
# Line 1678  class HomogeneousSaddlePointProblem(obje Line 1725  class HomogeneousSaddlePointProblem(obje
1725           :return: norm of v           :return: norm of v
1726           :rtype: non-negative ``float``           :rtype: non-negative ``float``
1727           """           """
1728           raise NotImplementedError,"no norm of v implemented."           raise NotImplementedError("no norm of v implemented.")
1729        def getDV(self, p, v, tol):        def getDV(self, p, v, tol):
1730           """           """
1731           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 1735  class HomogeneousSaddlePointProblem(obje
1735           :return: dv given as *dv= A^{-1} (f-A v-B^*p)*           :return: dv given as *dv= A^{-1} (f-A v-B^*p)*
1736           :note: Only *A* may depend on *v* and *p*           :note: Only *A* may depend on *v* and *p*
1737           """           """
1738           raise NotImplementedError,"no dv calculation implemented."           raise NotImplementedError("no dv calculation implemented.")
1739    
1740                    
1741        def Bv(self,v, tol):        def Bv(self,v, tol):
# Line 1698  class HomogeneousSaddlePointProblem(obje Line 1745  class HomogeneousSaddlePointProblem(obje
1745          :rtype: equal to the type of p          :rtype: equal to the type of p
1746          :note: boundary conditions on p should be zero!          :note: boundary conditions on p should be zero!
1747          """          """
1748          raise NotImplementedError, "no operator B implemented."          raise NotImplementedError("no operator B implemented.")
1749    
1750        def norm_Bv(self,Bv):        def norm_Bv(self,Bv):
1751          """          """
# Line 1707  class HomogeneousSaddlePointProblem(obje Line 1754  class HomogeneousSaddlePointProblem(obje
1754          :rtype: equal to the type of p          :rtype: equal to the type of p
1755          :note: boundary conditions on p should be zero!          :note: boundary conditions on p should be zero!
1756          """          """
1757          raise NotImplementedError, "no norm of Bv implemented."          raise NotImplementedError("no norm of Bv implemented.")
1758    
1759        def solve_AinvBt(self,dp, tol):        def solve_AinvBt(self,dp, tol):
1760           """           """
# Line 1717  class HomogeneousSaddlePointProblem(obje Line 1764  class HomogeneousSaddlePointProblem(obje
1764           :return: the solution of *A dv=B^*dp*           :return: the solution of *A dv=B^*dp*
1765           :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.
1766           """           """
1767           raise NotImplementedError,"no operator A implemented."           raise NotImplementedError("no operator A implemented.")
1768    
1769        def solve_prec(self,Bv, tol):        def solve_prec(self,Bv, tol):
1770           """           """
# Line 1726  class HomogeneousSaddlePointProblem(obje Line 1773  class HomogeneousSaddlePointProblem(obje
1773           :rtype: equal to the type of p           :rtype: equal to the type of p
1774           :note: boundary conditions on p should be zero!           :note: boundary conditions on p should be zero!
1775           """           """
1776           raise NotImplementedError,"no preconditioner for Schur complement implemented."           raise NotImplementedError("no preconditioner for Schur complement implemented.")
1777        #=============================================================        #=============================================================
1778        def __Aprod_PCG(self,dp):        def __Aprod_PCG(self,dp):
1779            dv=self.solve_AinvBt(dp, self.__subtol)            dv=self.solve_AinvBt(dp, self.__subtol)
# Line 1753  class HomogeneousSaddlePointProblem(obje Line 1800  class HomogeneousSaddlePointProblem(obje
1800            :rtype: ``float``            :rtype: ``float``
1801            """            """
1802            f=self.inner_p(p,p)            f=self.inner_p(p,p)
1803            if f<0: raise ValueError,"negative pressure norm."            if f<0: raise ValueError("negative pressure norm.")
1804            return math.sqrt(f)            return math.sqrt(f)
1805                
1806        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 1834  class HomogeneousSaddlePointProblem(obje
1834           self.verbose=verbose           self.verbose=verbose
1835           rtol=self.getTolerance()           rtol=self.getTolerance()
1836           atol=self.getAbsoluteTolerance()           atol=self.getAbsoluteTolerance()
1837    
1838             K_p=self.__K_p
1839             K_v=self.__K_v
1840           correction_step=0           correction_step=0
1841           converged=False           converged=False
          error=None  
1842           chi=None           chi=None
1843           gamma=self.__gamma           eps=None
1844           C_p=self.__C_p  
1845           C_v=self.__C_v           if self.verbose: print(("HomogeneousSaddlePointProblem: start iteration: rtol= %e, atol=%e"%(rtol, atol)))
1846           while not converged:           while not converged:
1847                if error== None or chi == None:  
1848                    gamma_new=gamma/self.__omega_conv               # get tolerance for velecity increment:
1849                else:               if chi == None:
1850                   if chi < self.__chi_max:                  rtol_v=self.__rtol_max
1851                      gamma_new=min(max(gamma*self.__omega_conv,1-chi*error/(self.__safety_factor*ATOL)), 1-chi/self.__chi_max)               else:
1852                   else:                  rtol_v=min(chi/K_v,self.__rtol_max)
1853                      gamma_new=gamma*self.__omega_div               rtol_v=max(rtol_v, self.__rtol_min)
1854                gamma=max(gamma_new, self.__gamma_min)               if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: rtol_v= %e"%(correction_step,rtol_v)))
1855                # calculate velocity for current pressure:               # get velocity increment:
1856                rtol_v=min(max(gamma/(1.+gamma)/C_v,self.__rtol_min), self.__rtol_max)               dv1=self.getDV(p,v,rtol_v)
1857                rtol_p=min(max(gamma/C_p, self.__rtol_min), self.__rtol_max)               v1=v+dv1
1858                self.__subtol=rtol_p**2               Bv1=self.Bv(v1, rtol_v)
1859                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)
1860                if self.verbose: print "HomogeneousSaddlePointProblem: subtolerance: %e"%self.__subtol               norm_dv1=self.norm_v(dv1)
1861                # 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)))
1862                dv1=self.getDV(p,v,rtol_v)               if norm_dv1*self.__theta < norm_Bv1:
1863                v1=v+dv1                  # get tolerance for pressure increment:
1864                Bv1=self.Bv(v1, self.__subtol)                  large_Bv1=True
1865                norm_Bv1=self.norm_Bv(Bv1)                  if chi == None or eps == None:
1866                norm_dv1=self.norm_v(dv1)                     rtol_p=self.__rtol_max
1867                norm_v1=self.norm_v(v1)                  else:
1868                ATOL=norm_v1*rtol+atol                     rtol_p=min(chi**2*eps/K_p/norm_Bv1, self.__rtol_max)
1869                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)
1870                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)))
1871                if max(norm_Bv1, norm_dv1) <= ATOL:                  # now we solve for the pressure increment dp from B*A^{-1}B^* dp = Bv1
1872                    converged=True                  if usePCG:
                   v=v1  
               else:  
                   # now we solve for the pressure increment dp from B*A^{-1}B^* dp = Bv1  
                   if usePCG:  
1873                      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)
1874                      v2=r[0]                      v2=r[0]
1875                      Bv2=r[1]                      Bv2=r[1]
1876                    else:                  else:
1877                        # don't use!!!!
1878                      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)
1879                      dv2=self.solve_AinvBt(dp, self.__subtol)                      dv2=self.solve_AinvBt(dp, self.__subtol)
1880                      v2=v1-dv2                      v2=v1-dv2
1881                      Bv2=self.Bv(v2, self.__subtol)                      Bv2=self.Bv(v2, self.__subtol)
1882                    #                  p2=p+dp
1883                    # convergence indicators:               else:
1884                    #                  large_Bv1=False
1885                    norm_v2=self.norm_v(v2)                  v2=v1
1886                    norm_dv2=self.norm_v(v2-v)                  p2=p
1887                    error_new=max(norm_dv2, norm_Bv1)               # update business:
1888                    norm_Bv2=self.norm_Bv(Bv2)               norm_dv2=self.norm_v(v2-v)
1889                    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)
1890                    if error !=None:               if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: v2 = %e, norm_dv2 = %e"%(correction_step, norm_v2, self.norm_v(v2-v))))
1891                        chi_new=error_new/error               eps, eps_old = max(norm_Bv1, norm_dv2), eps
1892                        if self.verbose: print "HomogeneousSaddlePointProblem: step %s: convergence rate = %e, est. error = %e"%(correction_step,chi_new, error_new)               if eps_old == None:
1893                        if chi != None:                    chi, chi_old = None, chi
1894                            gamma0=max(gamma, 1-chi/chi_new)               else:
1895                            C_p*=gamma0/gamma                    chi, chi_old = min(eps/ eps_old, self.__chi_max), chi
1896                            C_v*=gamma0/gamma*(1+gamma)/(1+gamma0)               if eps != None:
1897                        chi=chi_new                   if chi !=None:
1898                    else:                      if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: convergence rate = %e, correction = %e"%(correction_step,chi, eps)))
1899                        if self.verbose: print "HomogeneousSaddlePointProblem: step %s: est. error = %e"%(correction_step, error_new)                   else:
1900                        if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: correction = %e"%(correction_step, eps)))
1901                    error = error_new               if eps <= rtol*norm_v2+atol :
1902                    correction_step+=1                   converged = True
1903                    if correction_step>max_correction_steps:               else:
1904                          raise CorrectionFailed,"Given up after %d correction steps."%correction_step                   if correction_step>=max_correction_steps:
1905                    v,p=v2,p+dp                        raise CorrectionFailed("Given up after %d correction steps."%correction_step)
1906           if self.verbose: print "HomogeneousSaddlePointProblem: tolerance reached after %s steps."%correction_step                   if chi_old!=None:
1907       return v,p                      K_p=max(1,self.__reduction_factor*K_p,(chi-chi_old)/chi_old**2*K_p)
1908                        K_v=max(1,self.__reduction_factor*K_v,(chi-chi_old)/chi_old**2*K_p)
1909                        if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: new adjustment factor K = %e"%(correction_step,K_p)))
1910                 correction_step+=1
1911                 v,p =v2, p2
1912             if self.verbose: print(("HomogeneousSaddlePointProblem: tolerance reached after %s steps."%correction_step))
1913             return v,p
1914        #========================================================================        #========================================================================
1915        def setTolerance(self,tolerance=1.e-4):        def setTolerance(self,tolerance=1.e-4):
1916           """           """
# Line 1869  class HomogeneousSaddlePointProblem(obje Line 1920  class HomogeneousSaddlePointProblem(obje
1920           :type tolerance: non-negative ``float``           :type tolerance: non-negative ``float``
1921           """           """
1922           if tolerance<0:           if tolerance<0:
1923               raise ValueError,"tolerance must be positive."               raise ValueError("tolerance must be positive.")
1924           self.__rtol=tolerance           self.__rtol=tolerance
1925    
1926        def getTolerance(self):        def getTolerance(self):
# Line 1889  class HomogeneousSaddlePointProblem(obje Line 1940  class HomogeneousSaddlePointProblem(obje
1940           :type tolerance: non-negative ``float``           :type tolerance: non-negative ``float``
1941           """           """
1942           if tolerance<0:           if tolerance<0:
1943               raise ValueError,"tolerance must be non-negative."               raise ValueError("tolerance must be non-negative.")
1944           self.__atol=tolerance           self.__atol=tolerance
1945    
1946        def getAbsoluteTolerance(self):        def getAbsoluteTolerance(self):
# Line 1918  def MaskFromBoundaryTag(domain,*tags): Line 1969  def MaskFromBoundaryTag(domain,*tags):
1969     :rtype: `escript.Data` of rank 0     :rtype: `escript.Data` of rank 0
1970     """     """
1971     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1972     d=escript.Scalar(0.,escript.FunctionOnBoundary(domain))     d=escore.Scalar(0.,escore.FunctionOnBoundary(domain))
1973     for t in tags: d.setTaggedValue(t,1.)     for t in tags: d.setTaggedValue(t,1.)
1974     pde.setValue(y=d)     pde.setValue(y=d)
1975     return util.whereNonZero(pde.getRightHandSide())     return util.whereNonZero(pde.getRightHandSide())
# Line 1939  def MaskFromTag(domain,*tags): Line 1990  def MaskFromTag(domain,*tags):
1990     :rtype: `escript.Data` of rank 0     :rtype: `escript.Data` of rank 0
1991     """     """
1992     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1993     d=escript.Scalar(0.,escript.Function(domain))     d=escore.Scalar(0.,escore.Function(domain))
1994     for t in tags: d.setTaggedValue(t,1.)     for t in tags: d.setTaggedValue(t,1.)
1995     pde.setValue(Y=d)     pde.setValue(Y=d)
1996     return util.whereNonZero(pde.getRightHandSide())     return util.whereNonZero(pde.getRightHandSide())

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

  ViewVC Help
Powered by ViewVC 1.1.26