/[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 4446 by caltinay, Tue Jun 11 04:00:15 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 escriptcpp as escore
44  import linearPDEs  from . import linearPDEs
45  import numpy  from . import util
 import util  
46  import math  import math
47    import numpy
48    
49  class TimeIntegrationManager:  class TimeIntegrationManager(object):
50    """    """
51    A simple mechanism to manage time dependend values.    A simple mechanism to manage time dependend values.
52    
# 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 125  class TimeIntegrationManager: Line 126  class TimeIntegrationManager:
126           return out           return out
127    
128    
129  class Projector:  class Projector(object):
130    """    """
131    The Projector is a factory which projects a discontinuous function onto a    The Projector is a factory which projects a discontinuous function onto a
132    continuous function on a given domain.    continuous function on a given domain.
# 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.
168    
169      :param input_data: the data to be projected      :param input_data: the data to be projected
170      """      """
171      out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())      out=escore.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())
172      self.__pde.setValue(Y = escript.Data(), Y_reduced = escript.Data())      self.__pde.setValue(Y = escore.Data(), Y_reduced = escore.Data())
173      if input_data.getRank()==0:      if input_data.getRank()==0:
174          self.__pde.setValue(Y = input_data)          self.__pde.setValue(Y = input_data)
175          out=self.__pde.getSolution()          out=self.__pde.getSolution()
# Line 188  class Projector: Line 197  class Projector:
197                      out[i0,i1,i2,i3]=self.__pde.getSolution()                      out[i0,i1,i2,i3]=self.__pde.getSolution()
198      return out      return out
199    
200  class NoPDE:  class NoPDE(object):
201       """       """
202       Solves the following problem for u:       Solves the following problem for u:
203    
# Line 239  class NoPDE: Line 248  class NoPDE:
248           self.__q=q           self.__q=q
249           self.__r=r           self.__r=r
250           self.__u=None           self.__u=None
251           self.__function_space=escript.Solution(self.__domain)           self.__function_space=escore.Solution(self.__domain)
252    
253       def setReducedOn(self):       def setReducedOn(self):
254           """           """
255           Sets the `FunctionSpace` of the solution to `ReducedSolution`.           Sets the `FunctionSpace` of the solution to `ReducedSolution`.
256           """           """
257           self.__function_space=escript.ReducedSolution(self.__domain)           self.__function_space=escore.ReducedSolution(self.__domain)
258           self.__u=None           self.__u=None
259    
260       def setReducedOff(self):       def setReducedOff(self):
261           """           """
262           Sets the `FunctionSpace` of the solution to `Solution`.           Sets the `FunctionSpace` of the solution to `Solution`.
263           """           """
264           self.__function_space=escript.Solution(self.__domain)           self.__function_space=escore.Solution(self.__domain)
265           self.__u=None           self.__u=None
266    
267       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 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=escore.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=escore.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(escore.Data(self.__q,self.__function_space))
313                  self.__u*=(1.-q)                  self.__u*=(1.-q)
314                  if not self.__r==None: self.__u+=q*self.__r                  if not self.__r==None: self.__u+=q*self.__r
315           return self.__u           return self.__u
316    
317  class Locator:  class Locator(object):
318       """       """
319       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
320       coordinate x.       coordinate x.
# Line 325  class Locator: Line 334  class Locator:
334         :param x: location(s) of the Locator         :param x: location(s) of the Locator
335         :type x: ``numpy.ndarray`` or ``list`` of ``numpy.ndarray``         :type x: ``numpy.ndarray`` or ``list`` of ``numpy.ndarray``
336         """         """
337         if isinstance(where,escript.FunctionSpace):         if isinstance(where,escore.FunctionSpace):
338            self.__function_space=where            self.__function_space=where
339         else:         else:
340            self.__function_space=escript.ContinuousFunction(where)            self.__function_space=escore.ContinuousFunction(where)
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 400  class Locator: Line 410  class Locator:
410          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`
411          object otherwise the object is returned.          object otherwise the object is returned.
412          """          """
413          if isinstance(data,escript.Data):          if isinstance(data,escore.Data):
414             dat=util.interpolate(data,self.getFunctionSpace())             dat=util.interpolate(data,self.getFunctionSpace())
415             id=self.getId()             id=self.getId()
416             r=data.getRank()             r=data.getRank()
# 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, escore.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):
454      """      """
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, escore.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()
462      x_min=x.getTupleForGlobalDataPoint(*loc)      x_min=x.getTupleForGlobalDataPoint(*loc)
463      return Locator(arg.getFunctionSpace(),x_min)      return Locator(arg.getFunctionSpace(),x_min)
# Line 439  def getSupLocator(arg): Line 466  def getSupLocator(arg):
466      """      """
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, escore.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()
474      x_min=x.getTupleForGlobalDataPoint(*loc)      x_min=x.getTupleForGlobalDataPoint(*loc)
475      return Locator(arg.getFunctionSpace(),x_min)      return Locator(arg.getFunctionSpace(),x_min)
476                
477    
478  class SolverSchemeException(Exception):  class SolverSchemeException(Exception):
479     """     """
# 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_
748              #              #
         #   adjust subtol_  
         #  
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
758              #     if iter_restart -1 is returned as sub_iter              #     if iter_restart -1 is returned as sub_iter
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 1316  class ArithmeticTuple(object): Line 1343  class ArithmeticTuple(object):
1343    
1344         from esys.escript import Data         from esys.escript import Data
1345         from numpy import array         from numpy import array
1346         a=Data(...)         a=eData(...)
1347         b=array([1.,4.])         b=array([1.,4.])
1348         x=ArithmeticTuple(a,b)         x=ArithmeticTuple(a,b)
1349         y=5.*x         y=5.*x
# 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):
1395                    if self.__isEmpty(self[i]) or self.__isEmpty(other[i]):
1396                        out.append(escore.Data())
1397                    else:
1398                        out.append(self[i]*other[i])
1399         except TypeError:         except TypeError:
1400             for i in range(len(self)): out.append(self[i]*other)             for i in range(len(self)):  
1401                    if self.__isEmpty(self[i]) or self.__isEmpty(other):
1402                        out.append(escore.Data())
1403                    else:
1404                        out.append(self[i]*other)
1405         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1406    
1407     def __rmul__(self,other):     def __rmul__(self,other):
1408         """        """
1409         Scales by ``other`` from the left.        Scales by ``other`` from the left.
1410    
1411         :param other: scaling factor        :param other: scaling factor
1412         :type other: ``float``        :type other: ``float``
1413         :return: itemwise other*self        :return: itemwise other*self
1414         :rtype: `ArithmeticTuple`        :rtype: `ArithmeticTuple`
1415         """        """
1416         out=[]        out=[]
1417         try:        try:
1418             l=len(other)            l=len(other)
1419             if l!=len(self):            if l!=len(self):
1420                 raise ValueError,"length of arguments don't match."                raise ValueError("length of arguments don't match.")
1421             for i in range(l): out.append(other[i]*self[i])            for i in range(l):
1422         except TypeError:                  if self.__isEmpty(self[i]) or self.__isEmpty(other[i]):
1423             for i in range(len(self)): out.append(other*self[i])                      out.append(escore.Data())
1424         return ArithmeticTuple(*tuple(out))                  else:
1425                        out.append(other[i]*self[i])
1426          except TypeError:
1427              for i in range(len(self)):  
1428                    if self.__isEmpty(self[i]) or self.__isEmpty(other):
1429                        out.append(escore.Data())
1430                    else:
1431                        out.append(other*self[i])
1432          return ArithmeticTuple(*tuple(out))
1433    
1434     def __div__(self,other):     def __div__(self,other):
1435         """         """
# Line 1400  class ArithmeticTuple(object): Line 1443  class ArithmeticTuple(object):
1443         return self*(1/other)         return self*(1/other)
1444    
1445     def __rdiv__(self,other):     def __rdiv__(self,other):
1446         """        """
1447         Scales by (1/``other``) from the left.        Scales by (1/``other``) from the left.
1448    
1449         :param other: scaling factor        :param other: scaling factor
1450         :type other: ``float``        :type other: ``float``
1451         :return: itemwise other/self        :return: itemwise other/self
1452         :rtype: `ArithmeticTuple`        :rtype: `ArithmeticTuple`
1453         """        """
1454         out=[]        out=[]
1455         try:        try:
1456             l=len(other)            l=len(other)
1457             if l!=len(self):            if l!=len(self):
1458                 raise ValueError,"length of arguments don't match."                raise ValueError("length of arguments don't match.")
1459             for i in range(l): out.append(other[i]/self[i])            
1460         except TypeError:            for i in range(l):
1461             for i in range(len(self)): out.append(other/self[i])                  if self.__isEmpty(self[i]):
1462         return ArithmeticTuple(*tuple(out))                      raise ZeroDivisionError("in component %s"%i)
1463                    else:
1464                        if self.__isEmpty(other[i]):
1465                            out.append(escore.Data())
1466                        else:
1467                            out.append(other[i]/self[i])
1468          except TypeError:
1469              for i in range(len(self)):
1470                    if self.__isEmpty(self[i]):
1471                        raise ZeroDivisionError("in component %s"%i)
1472                    else:
1473                        if self.__isEmpty(other):
1474                            out.append(escore.Data())
1475                        else:
1476                            out.append(other/self[i])
1477          return ArithmeticTuple(*tuple(out))
1478    
1479     def __iadd__(self,other):     def __iadd__(self,other):
1480         """        """
1481         Inplace addition of ``other`` to self.        Inplace addition of ``other`` to self.
1482    
1483         :param other: increment        :param other: increment
1484         :type other: ``ArithmeticTuple``        :type other: ``ArithmeticTuple``
1485         """        """
1486         if len(self) != len(other):        if len(self) != len(other):
1487             raise ValueError,"tuple lengths must match."            raise ValueError("tuple lengths must match.")
1488         for i in range(len(self)):        for i in range(len(self)):
1489             self.__items[i]+=other[i]            if self.__isEmpty(self.__items[i]):
1490         return self                self.__items[i]=other[i]
1491              else:
1492                  self.__items[i]+=other[i]
1493                  
1494          return self
1495    
1496     def __add__(self,other):     def __add__(self,other):
1497         """        """
1498         Adds ``other`` to self.        Adds ``other`` to self.
1499    
1500         :param other: increment        :param other: increment
1501         :type other: ``ArithmeticTuple``        :type other: ``ArithmeticTuple``
1502         """        """
1503         out=[]        out=[]
1504         try:        try:
1505             l=len(other)            l=len(other)
1506             if l!=len(self):            if l!=len(self):
1507                 raise ValueError,"length of arguments don't match."                raise ValueError("length of arguments don't match.")
1508             for i in range(l): out.append(self[i]+other[i])            for i in range(l):
1509         except TypeError:                  if self.__isEmpty(self[i]):
1510             for i in range(len(self)): out.append(self[i]+other)                      out.append(other[i])
1511         return ArithmeticTuple(*tuple(out))                  elif self.__isEmpty(other[i]):
1512                        out.append(self[i])
1513                    else:
1514                        out.append(self[i]+other[i])
1515          except TypeError:
1516                for i in range(len(self)):    
1517                    if self.__isEmpty(self[i]):
1518                        out.append(other)
1519                    elif self.__isEmpty(other):
1520                        out.append(self[i])
1521                    else:
1522                        out.append(self[i]+other)
1523          return ArithmeticTuple(*tuple(out))
1524    
1525     def __sub__(self,other):     def __sub__(self,other):
1526         """        """
1527         Subtracts ``other`` from self.        Subtracts ``other`` from self.
1528    
1529         :param other: decrement        :param other: decrement
1530         :type other: ``ArithmeticTuple``        :type other: ``ArithmeticTuple``
1531         """        """
1532         out=[]        out=[]
1533         try:        try:
1534             l=len(other)            l=len(other)
1535             if l!=len(self):            if l!=len(self):
1536                 raise ValueError,"length of arguments don't match."                raise ValueError("length of arguments don't match.")
1537             for i in range(l): out.append(self[i]-other[i])            for i in range(l):
1538         except TypeError:                  if self.__isEmpty(other[i]):
1539             for i in range(len(self)): out.append(self[i]-other)                      out.append(self[i])
1540         return ArithmeticTuple(*tuple(out))                  elif self.__isEmpty(self[i]):
1541                        out.append(-other[i])
1542                    else:
1543                        out.append(self[i]-other[i])
1544          except TypeError:
1545                for i in range(len(self)):    
1546                    if  self.__isEmpty(other):
1547                        out.append(self[i])
1548                    elif self.__isEmpty(self[i]):
1549                        out.append(-other)
1550                    else:
1551                        out.append(self[i]-other)
1552                        
1553          return ArithmeticTuple(*tuple(out))
1554    
1555     def __isub__(self,other):     def __isub__(self,other):
1556         """        """
1557         Inplace subtraction of ``other`` from self.        Inplace subtraction of ``other`` from self.
1558    
1559         :param other: decrement        :param other: decrement
1560         :type other: ``ArithmeticTuple``        :type other: ``ArithmeticTuple``
1561         """        """
1562         if len(self) != len(other):        if len(self) != len(other):
1563             raise ValueError,"tuple length must match."            raise ValueError("tuple length must match.")
1564         for i in range(len(self)):        for i in range(len(self)):
1565             self.__items[i]-=other[i]            if not self.__isEmpty(other[i]):
1566         return self                if self.__isEmpty(self.__items[i]):
1567                      self.__items[i]=-other[i]
1568                  else:
1569                      self.__items[i]=other[i]
1570          return self
1571    
1572     def __neg__(self):     def __neg__(self):
1573         """        """
1574         Negates values.        Negates values.
1575         """        """
1576         out=[]        out=[]
1577         for i in range(len(self)):        for i in range(len(self)):
1578             out.append(-self[i])            if self.__isEmpty(self[i]):
1579         return ArithmeticTuple(*tuple(out))                out.append(escore.Data())
1580              else:
1581                  out.append(-self[i])
1582              
1583          return ArithmeticTuple(*tuple(out))
1584       def __isEmpty(self, d):
1585        if isinstance(d, escore.Data):
1586            return d.isEmpty()
1587        else:
1588            return False
1589                    
1590       def __str__(self):
1591        s="("
1592        for i in self:
1593          s=s+str(i)+", "
1594        s=s+")"
1595        return s
1596    
1597  class HomogeneousSaddlePointProblem(object):  class HomogeneousSaddlePointProblem(object):
1598        """        """
# Line 1501  class HomogeneousSaddlePointProblem(obje Line 1607  class HomogeneousSaddlePointProblem(obje
1607        *A* may depend weakly on *v* and *p*.        *A* may depend weakly on *v* and *p*.
1608        """        """
1609        def __init__(self, **kwargs):        def __init__(self, **kwargs):
1610      """          """
1611      initializes the saddle point problem          initializes the saddle point problem
1612      """          """
1613          self.resetControlParameters()          self.resetControlParameters()
1614          self.setTolerance()          self.setTolerance()
1615          self.setAbsoluteTolerance()          self.setAbsoluteTolerance()
1616        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):
1617           """           """
1618           sets a control parameter           sets a control parameter
1619    
1620           :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
1621           :type gamma: ``float``           :type K_p: ``float``
1622           :param gamma_min: minimum value for ``gamma``.           :param K_v: initial value for constant to adjust velocity tolerance
1623           :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``  
1624           :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.
1625           :type rtol_max: ``float``           :type rtol_max: ``float``
1626           :param chi: initial convergence measure.           :param chi_max: maximum tolerable converegence rate.
1627           :type chi: ``float``           :type chi_max: ``float``
1628           :param C_p: initial value for constant to adjust pressure tolerance           :param reduction_factor: reduction factor for adjustment factors.
1629           :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``  
1630           """           """
1631           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)
1632    
1633        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):
1634           """           """
1635           sets a control parameter           sets a control parameter
1636    
1637           :param gamma: ``1/(1-gamma)`` controls the perturbation of the converegence rate due to termination errors in the subproblems.  
1638           :type gamma: ``float``           :param K_p: initial value for constant to adjust pressure tolerance
1639           :param gamma_min: minimum value for ``gamma``.           :type K_p: ``float``
1640           :type gamma_min: ``float``           :param K_v: initial value for constant to adjust velocity tolerance
1641           :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``  
1642           :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.
1643           :type rtol_max: ``float``           :type rtol_max: ``float``
1644           :param chi: initial convergence measure.           :param chi_max: maximum tolerable converegence rate.
1645           :type chi: ``float``           :type chi_max: ``float``
1646           :param C_p: initial value for constant to adjust pressure tolerance           :type reduction_factor: ``float``
1647           :type C_p: ``float``           """
1648           :param C_v: initial value for constant to adjust velocity tolerance           if not K_p == None:
1649           :type C_v: ``float``              if K_p<1:
1650           :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."  
1651           else:           else:
1652              omega_div=self.__omega_div              K_p=self.__K_p
1653    
1654           if not omega_conv == None:           if not K_v == None:
1655              if omega_conv<1:              if K_v<1:
1656                 raise ValueError,"omega_conv needs to be greater or equal to 1."                 raise ValueError("K_v need to be greater or equal to 1.")
1657           else:           else:
1658              omega_conv=self.__omega_conv              K_v=self.__K_v
   
          if not rtol_min == None:  
             if rtol_min<=0 or rtol_min>=1:  
                raise ValueError,"rtol_min needs to be positive and less than 1."  
          else:  
             rtol_min=self.__rtol_min  
1659    
1660           if not rtol_max == None:           if not rtol_max == None:
1661              if rtol_max<=0 or rtol_max>=1:              if rtol_max<=0 or rtol_max>=1:
1662                 raise ValueError,"rtol_max needs to be positive and less than 1."                 raise ValueError("rtol_max needs to be positive and less than 1.")
1663           else:           else:
1664              rtol_max=self.__rtol_max              rtol_max=self.__rtol_max
1665    
1666           if not chi == None:           if not rtol_min == None:
1667              if chi<=0 or chi>1:              if rtol_min<=0 or rtol_min>=1:
1668                 raise ValueError,"chi needs to be positive and less than 1."                 raise ValueError("rtol_min needs to be positive and less than 1.")
1669           else:           else:
1670              chi=self.__chi              rtol_min=self.__rtol_min
1671    
1672           if not C_p == None:           if not chi_max == None:
1673              if C_p<1:              if chi_max<=0 or chi_max>=1:
1674                 raise ValueError,"C_p need to be greater or equal to 1."                 raise ValueError("chi_max needs to be positive and less than 1.")
1675           else:           else:
1676              C_p=self.__C_p              chi_max = self.__chi_max
1677    
1678           if not C_v == None:           if not reduction_factor == None:
1679              if C_v<1:              if reduction_factor<=0 or reduction_factor>1:
1680                 raise ValueError,"C_v need to be greater or equal to 1."                 raise ValueError("reduction_factor need to be between zero and one.")
1681           else:           else:
1682              C_v=self.__C_v              reduction_factor=self.__reduction_factor
1683    
1684           if not safety_factor == None:           if not theta == None:
1685              if safety_factor<=0 or safety_factor>1:              if theta<=0 or theta>1:
1686                 raise ValueError,"safety_factor need to be between zero and one."                 raise ValueError("theta need to be between zero and one.")
1687           else:           else:
1688              safety_factor=self.__safety_factor              theta=self.__theta
1689    
1690           if gamma<gamma_min:           if rtol_min>=rtol_max:
1691                 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  
1692           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  
1693           self.__rtol_max = rtol_max           self.__rtol_max = rtol_max
1694           self.__chi = chi           self.__K_p = K_p
1695           self.__C_p = C_p           self.__K_v = K_v
1696           self.__C_v = C_v           self.__reduction_factor = reduction_factor
1697           self.__safety_factor = safety_factor           self.__theta = theta
1698             self.__rtol_min=rtol_min
1699    
1700        #=============================================================        #=============================================================
1701        def inner_pBv(self,p,Bv):        def inner_pBv(self,p,Bv):
# Line 1657  class HomogeneousSaddlePointProblem(obje Line 1708  class HomogeneousSaddlePointProblem(obje
1708           :rtype: ``float``           :rtype: ``float``
1709           :note: used if PCG is applied.           :note: used if PCG is applied.
1710           """           """
1711           raise NotImplementedError,"no inner product for p and Bv implemented."           raise NotImplementedError("no inner product for p and Bv implemented.")
1712    
1713        def inner_p(self,p0,p1):        def inner_p(self,p0,p1):
1714           """           """
# Line 1668  class HomogeneousSaddlePointProblem(obje Line 1719  class HomogeneousSaddlePointProblem(obje
1719           :return: inner product of p0 and p1           :return: inner product of p0 and p1
1720           :rtype: ``float``           :rtype: ``float``
1721           """           """
1722           raise NotImplementedError,"no inner product for p implemented."           raise NotImplementedError("no inner product for p implemented.")
1723        
1724        def norm_v(self,v):        def norm_v(self,v):
1725           """           """
# Line 1678  class HomogeneousSaddlePointProblem(obje Line 1729  class HomogeneousSaddlePointProblem(obje
1729           :return: norm of v           :return: norm of v
1730           :rtype: non-negative ``float``           :rtype: non-negative ``float``
1731           """           """
1732           raise NotImplementedError,"no norm of v implemented."           raise NotImplementedError("no norm of v implemented.")
1733        def getDV(self, p, v, tol):        def getDV(self, p, v, tol):
1734           """           """
1735           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 1739  class HomogeneousSaddlePointProblem(obje
1739           :return: dv given as *dv= A^{-1} (f-A v-B^*p)*           :return: dv given as *dv= A^{-1} (f-A v-B^*p)*
1740           :note: Only *A* may depend on *v* and *p*           :note: Only *A* may depend on *v* and *p*
1741           """           """
1742           raise NotImplementedError,"no dv calculation implemented."           raise NotImplementedError("no dv calculation implemented.")
1743    
1744                    
1745        def Bv(self,v, tol):        def Bv(self,v, tol):
# Line 1698  class HomogeneousSaddlePointProblem(obje Line 1749  class HomogeneousSaddlePointProblem(obje
1749          :rtype: equal to the type of p          :rtype: equal to the type of p
1750          :note: boundary conditions on p should be zero!          :note: boundary conditions on p should be zero!
1751          """          """
1752          raise NotImplementedError, "no operator B implemented."          raise NotImplementedError("no operator B implemented.")
1753    
1754        def norm_Bv(self,Bv):        def norm_Bv(self,Bv):
1755          """          """
# Line 1707  class HomogeneousSaddlePointProblem(obje Line 1758  class HomogeneousSaddlePointProblem(obje
1758          :rtype: equal to the type of p          :rtype: equal to the type of p
1759          :note: boundary conditions on p should be zero!          :note: boundary conditions on p should be zero!
1760          """          """
1761          raise NotImplementedError, "no norm of Bv implemented."          raise NotImplementedError("no norm of Bv implemented.")
1762    
1763        def solve_AinvBt(self,dp, tol):        def solve_AinvBt(self,dp, tol):
1764           """           """
# Line 1717  class HomogeneousSaddlePointProblem(obje Line 1768  class HomogeneousSaddlePointProblem(obje
1768           :return: the solution of *A dv=B^*dp*           :return: the solution of *A dv=B^*dp*
1769           :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.
1770           """           """
1771           raise NotImplementedError,"no operator A implemented."           raise NotImplementedError("no operator A implemented.")
1772    
1773        def solve_prec(self,Bv, tol):        def solve_prec(self,Bv, tol):
1774           """           """
# Line 1726  class HomogeneousSaddlePointProblem(obje Line 1777  class HomogeneousSaddlePointProblem(obje
1777           :rtype: equal to the type of p           :rtype: equal to the type of p
1778           :note: boundary conditions on p should be zero!           :note: boundary conditions on p should be zero!
1779           """           """
1780           raise NotImplementedError,"no preconditioner for Schur complement implemented."           raise NotImplementedError("no preconditioner for Schur complement implemented.")
1781        #=============================================================        #=============================================================
1782        def __Aprod_PCG(self,dp):        def __Aprod_PCG(self,dp):
1783            dv=self.solve_AinvBt(dp, self.__subtol)            dv=self.solve_AinvBt(dp, self.__subtol)
# Line 1753  class HomogeneousSaddlePointProblem(obje Line 1804  class HomogeneousSaddlePointProblem(obje
1804            :rtype: ``float``            :rtype: ``float``
1805            """            """
1806            f=self.inner_p(p,p)            f=self.inner_p(p,p)
1807            if f<0: raise ValueError,"negative pressure norm."            if f<0: raise ValueError("negative pressure norm.")
1808            return math.sqrt(f)            return math.sqrt(f)
1809                    
1810        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):
1811           """           """
1812           Solves the saddle point problem using initial guesses v and p.           Solves the saddle point problem using initial guesses v and p.
# Line 1787  class HomogeneousSaddlePointProblem(obje Line 1838  class HomogeneousSaddlePointProblem(obje
1838           self.verbose=verbose           self.verbose=verbose
1839           rtol=self.getTolerance()           rtol=self.getTolerance()
1840           atol=self.getAbsoluteTolerance()           atol=self.getAbsoluteTolerance()
1841    
1842             K_p=self.__K_p
1843             K_v=self.__K_v
1844           correction_step=0           correction_step=0
1845           converged=False           converged=False
          error=None  
1846           chi=None           chi=None
1847           gamma=self.__gamma           eps=None
1848           C_p=self.__C_p  
1849           C_v=self.__C_v           if self.verbose: print(("HomogeneousSaddlePointProblem: start iteration: rtol= %e, atol=%e"%(rtol, atol)))
1850           while not converged:           while not converged:
1851                if error== None or chi == None:  
1852                    gamma_new=gamma/self.__omega_conv               # get tolerance for velecity increment:
1853                else:               if chi == None:
1854                   if chi < self.__chi_max:                  rtol_v=self.__rtol_max
1855                      gamma_new=min(max(gamma*self.__omega_conv,1-chi*error/(self.__safety_factor*ATOL)), 1-chi/self.__chi_max)               else:
1856                   else:                  rtol_v=min(chi/K_v,self.__rtol_max)
1857                      gamma_new=gamma*self.__omega_div               rtol_v=max(rtol_v, self.__rtol_min)
1858                gamma=max(gamma_new, self.__gamma_min)               if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: rtol_v= %e"%(correction_step,rtol_v)))
1859                # calculate velocity for current pressure:               # get velocity increment:
1860                rtol_v=min(max(gamma/(1.+gamma)/C_v,self.__rtol_min), self.__rtol_max)               dv1=self.getDV(p,v,rtol_v)
1861                rtol_p=min(max(gamma/C_p, self.__rtol_min), self.__rtol_max)               v1=v+dv1
1862                self.__subtol=rtol_p**2               Bv1=self.Bv(v1, rtol_v)
1863                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)
1864                if self.verbose: print "HomogeneousSaddlePointProblem: subtolerance: %e"%self.__subtol               norm_dv1=self.norm_v(dv1)
1865                # 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)))
1866                dv1=self.getDV(p,v,rtol_v)               if norm_dv1*self.__theta < norm_Bv1:
1867                v1=v+dv1                  # get tolerance for pressure increment:
1868                Bv1=self.Bv(v1, self.__subtol)                  large_Bv1=True
1869                norm_Bv1=self.norm_Bv(Bv1)                  if chi == None or eps == None:
1870                norm_dv1=self.norm_v(dv1)                     rtol_p=self.__rtol_max
1871                norm_v1=self.norm_v(v1)                  else:
1872                ATOL=norm_v1*rtol+atol                     rtol_p=min(chi**2*eps/K_p/norm_Bv1, self.__rtol_max)
1873                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)
1874                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)))
1875                if max(norm_Bv1, norm_dv1) <= ATOL:                  # now we solve for the pressure increment dp from B*A^{-1}B^* dp = Bv1
1876                    converged=True                  if usePCG:
                   v=v1  
               else:  
                   # now we solve for the pressure increment dp from B*A^{-1}B^* dp = Bv1  
                   if usePCG:  
1877                      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)
1878                      v2=r[0]                      v2=r[0]
1879                      Bv2=r[1]                      Bv2=r[1]
1880                    else:                  else:
1881                        # don't use!!!!
1882                      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)
1883                      dv2=self.solve_AinvBt(dp, self.__subtol)                      dv2=self.solve_AinvBt(dp, self.__subtol)
1884                      v2=v1-dv2                      v2=v1-dv2
1885                      Bv2=self.Bv(v2, self.__subtol)                      Bv2=self.Bv(v2, self.__subtol)
1886                    #                  p2=p+dp
1887                    # convergence indicators:               else:
1888                    #                  large_Bv1=False
1889                    norm_v2=self.norm_v(v2)                  v2=v1
1890                    norm_dv2=self.norm_v(v2-v)                  p2=p
1891                    error_new=max(norm_dv2, norm_Bv1)               # update business:
1892                    norm_Bv2=self.norm_Bv(Bv2)               norm_dv2=self.norm_v(v2-v)
1893                    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)
1894                    if error !=None:               if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: v2 = %e, norm_dv2 = %e"%(correction_step, norm_v2, self.norm_v(v2-v))))
1895                        chi_new=error_new/error               eps, eps_old = max(norm_Bv1, norm_dv2), eps
1896                        if self.verbose: print "HomogeneousSaddlePointProblem: step %s: convergence rate = %e, est. error = %e"%(correction_step,chi_new, error_new)               if eps_old == None:
1897                        if chi != None:                    chi, chi_old = None, chi
1898                            gamma0=max(gamma, 1-chi/chi_new)               else:
1899                            C_p*=gamma0/gamma                    chi, chi_old = min(eps/ eps_old, self.__chi_max), chi
1900                            C_v*=gamma0/gamma*(1+gamma)/(1+gamma0)               if eps != None:
1901                        chi=chi_new                   if chi !=None:
1902                    else:                      if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: convergence rate = %e, correction = %e"%(correction_step,chi, eps)))
1903                        if self.verbose: print "HomogeneousSaddlePointProblem: step %s: est. error = %e"%(correction_step, error_new)                   else:
1904                        if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: correction = %e"%(correction_step, eps)))
1905                    error = error_new               if eps <= rtol*norm_v2+atol :
1906                    correction_step+=1                   converged = True
1907                    if correction_step>max_correction_steps:               else:
1908                          raise CorrectionFailed,"Given up after %d correction steps."%correction_step                   if correction_step>=max_correction_steps:
1909                    v,p=v2,p+dp                        raise CorrectionFailed("Given up after %d correction steps."%correction_step)
1910           if self.verbose: print "HomogeneousSaddlePointProblem: tolerance reached after %s steps."%correction_step                   if chi_old!=None:
1911       return v,p                      K_p=max(1,self.__reduction_factor*K_p,(chi-chi_old)/chi_old**2*K_p)
1912                        K_v=max(1,self.__reduction_factor*K_v,(chi-chi_old)/chi_old**2*K_p)
1913                        if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: new adjustment factor K = %e"%(correction_step,K_p)))
1914                 correction_step+=1
1915                 v,p =v2, p2
1916             if self.verbose: print(("HomogeneousSaddlePointProblem: tolerance reached after %s steps."%correction_step))
1917             return v,p
1918        #========================================================================        #========================================================================
1919        def setTolerance(self,tolerance=1.e-4):        def setTolerance(self,tolerance=1.e-4):
1920           """           """
# Line 1869  class HomogeneousSaddlePointProblem(obje Line 1924  class HomogeneousSaddlePointProblem(obje
1924           :type tolerance: non-negative ``float``           :type tolerance: non-negative ``float``
1925           """           """
1926           if tolerance<0:           if tolerance<0:
1927               raise ValueError,"tolerance must be positive."               raise ValueError("tolerance must be positive.")
1928           self.__rtol=tolerance           self.__rtol=tolerance
1929    
1930        def getTolerance(self):        def getTolerance(self):
# Line 1889  class HomogeneousSaddlePointProblem(obje Line 1944  class HomogeneousSaddlePointProblem(obje
1944           :type tolerance: non-negative ``float``           :type tolerance: non-negative ``float``
1945           """           """
1946           if tolerance<0:           if tolerance<0:
1947               raise ValueError,"tolerance must be non-negative."               raise ValueError("tolerance must be non-negative.")
1948           self.__atol=tolerance           self.__atol=tolerance
1949    
1950        def getAbsoluteTolerance(self):        def getAbsoluteTolerance(self):
# Line 1918  def MaskFromBoundaryTag(domain,*tags): Line 1973  def MaskFromBoundaryTag(domain,*tags):
1973     :rtype: `escript.Data` of rank 0     :rtype: `escript.Data` of rank 0
1974     """     """
1975     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1976     d=escript.Scalar(0.,escript.FunctionOnBoundary(domain))     d=escore.Scalar(0.,escore.FunctionOnBoundary(domain))
1977     for t in tags: d.setTaggedValue(t,1.)     for t in tags: d.setTaggedValue(t,1.)
1978     pde.setValue(y=d)     pde.setValue(y=d)
1979     return util.whereNonZero(pde.getRightHandSide())     return util.whereNonZero(pde.getRightHandSide())
# Line 1939  def MaskFromTag(domain,*tags): Line 1994  def MaskFromTag(domain,*tags):
1994     :rtype: `escript.Data` of rank 0     :rtype: `escript.Data` of rank 0
1995     """     """
1996     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1997     d=escript.Scalar(0.,escript.Function(domain))     d=escore.Scalar(0.,escore.Function(domain))
1998     for t in tags: d.setTaggedValue(t,1.)     for t in tags: d.setTaggedValue(t,1.)
1999     pde.setValue(Y=d)     pde.setValue(Y=d)
2000     return util.whereNonZero(pde.getRightHandSide())     return util.whereNonZero(pde.getRightHandSide())

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

  ViewVC Help
Powered by ViewVC 1.1.26