/[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 782 by bcumming, Tue Jul 18 00:47:47 2006 UTC revision 2100 by gross, Wed Nov 26 08:13:00 2008 UTC
# Line 1  Line 1 
1  # $Id$  
2    ########################################################
3    #
4    # Copyright (c) 2003-2008 by University of Queensland
5    # Earth Systems Science Computational Center (ESSCC)
6    # http://www.uq.edu.au/esscc
7    #
8    # Primary Business: Queensland, Australia
9    # Licensed under the Open Software License version 3.0
10    # http://www.opensource.org/licenses/osl-3.0.php
11    #
12    ########################################################
13    
14    __copyright__="""Copyright (c) 2003-2008 by University of Queensland
15    Earth Systems Science Computational Center (ESSCC)
16    http://www.uq.edu.au/esscc
17    Primary Business: Queensland, Australia"""
18    __license__="""Licensed under the Open Software License version 3.0
19    http://www.opensource.org/licenses/osl-3.0.php"""
20    __url__="http://www.uq.edu.au/esscc/escript-finley"
21    
22  """  """
23  Provides some tools related to PDEs.  Provides some tools related to PDEs.
# Line 7  Currently includes: Line 26  Currently includes:
26      - Projector - to project a discontinuous      - Projector - to project a discontinuous
27      - Locator - to trace values in data objects at a certain location      - Locator - to trace values in data objects at a certain location
28      - TimeIntegrationManager - to handel extraplotion in time      - TimeIntegrationManager - to handel extraplotion in time
29        - SaddlePointProblem - solver for Saddle point problems using the inexact uszawa scheme
30    
31  @var __author__: name of author  @var __author__: name of author
32  @var __copyright__: copyrights  @var __copyright__: copyrights
# Line 17  Currently includes: Line 37  Currently includes:
37  """  """
38    
39  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
 __copyright__="""  Copyright (c) 2006 by ACcESS MNRF  
                     http://www.access.edu.au  
                 Primary Business: Queensland, Australia"""  
 __license__="""Licensed under the Open Software License version 3.0  
              http://www.opensource.org/licenses/osl-3.0.php"""  
 __url__="http://www.iservo.edu.au/esys"  
 __version__="$Revision$"  
 __date__="$Date$"  
40    
41    
42  import escript  import escript
43  import linearPDEs  import linearPDEs
44  import numarray  import numarray
45  import util  import util
46    import math
47    
48    ##### Added by Artak
49    # from Numeric import zeros,Int,Float64
50    ###################################
51    
52    
53  class TimeIntegrationManager:  class TimeIntegrationManager:
54    """    """
# Line 131  class Projector: Line 149  class Projector:
149      self.__pde.setValue(D = 1.)      self.__pde.setValue(D = 1.)
150      return      return
151    
   def __del__(self):  
     return  
   
152    def __call__(self, input_data):    def __call__(self, input_data):
153      """      """
154      Projects input_data onto a continuous function      Projects input_data onto a continuous function
# Line 141  class Projector: Line 156  class Projector:
156      @param input_data: The input_data to be projected.      @param input_data: The input_data to be projected.
157      """      """
158      out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())      out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())
159        self.__pde.setValue(Y = escript.Data(), Y_reduced = escript.Data())
160      if input_data.getRank()==0:      if input_data.getRank()==0:
161          self.__pde.setValue(Y = input_data)          self.__pde.setValue(Y = input_data)
162          out=self.__pde.getSolution()          out=self.__pde.getSolution()
# Line 296  class Locator: Line 312  class Locator:
312         Initializes a Locator to access values in Data objects on the Doamin         Initializes a Locator to access values in Data objects on the Doamin
313         or FunctionSpace where for the sample point which         or FunctionSpace where for the sample point which
314         closest to the given point x.         closest to the given point x.
315    
316           @param where: function space
317           @type where: L{escript.FunctionSpace}
318           @param x: coefficient of the solution.
319           @type x: L{numarray.NumArray} or C{list} of L{numarray.NumArray}
320         """         """
321         if isinstance(where,escript.FunctionSpace):         if isinstance(where,escript.FunctionSpace):
322            self.__function_space=where            self.__function_space=where
323         else:         else:
324            self.__function_space=escript.ContinuousFunction(where)            self.__function_space=escript.ContinuousFunction(where)
325         self.__id=util.length(self.__function_space.getX()-x[:self.__function_space.getDim()]).mindp()         if isinstance(x, list):
326               self.__id=[]
327               for p in x:
328                  self.__id.append(util.length(self.__function_space.getX()-p[:self.__function_space.getDim()]).minGlobalDataPoint())
329           else:
330               self.__id=util.length(self.__function_space.getX()-x[:self.__function_space.getDim()]).minGlobalDataPoint()
331    
332       def __str__(self):       def __str__(self):
333         """         """
334         Returns the coordinates of the Locator as a string.         Returns the coordinates of the Locator as a string.
335         """         """
336         return "<Locator %s>"%str(self.getX())         x=self.getX()
337           if instance(x,list):
338              out="["
339              first=True
340              for xx in x:
341                if not first:
342                    out+=","
343                else:
344                    first=False
345                out+=str(xx)
346              out+="]>"
347           else:
348              out=str(x)
349           return out
350    
351         def getX(self):
352            """
353        Returns the exact coordinates of the Locator.
354        """
355            return self(self.getFunctionSpace().getX())
356    
357       def getFunctionSpace(self):       def getFunctionSpace(self):
358          """          """
# Line 315  class Locator: Line 360  class Locator:
360      """      """
361          return self.__function_space          return self.__function_space
362    
363       def getId(self):       def getId(self,item=None):
364          """          """
365      Returns the identifier of the location.      Returns the identifier of the location.
366      """      """
367          return self.__id          if item == None:
368               return self.__id
369            else:
370               if isinstance(self.__id,list):
371                  return self.__id[item]
372               else:
373                  return self.__id
374    
      def getX(self):  
         """  
     Returns the exact coordinates of the Locator.  
     """  
         return self(self.getFunctionSpace().getX())  
375    
376       def __call__(self,data):       def __call__(self,data):
377          """          """
# Line 341  class Locator: Line 387  class Locator:
387      """      """
388          if isinstance(data,escript.Data):          if isinstance(data,escript.Data):
389             if data.getFunctionSpace()==self.getFunctionSpace():             if data.getFunctionSpace()==self.getFunctionSpace():
390               out=data.convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])               dat=data
              #out=data.convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1],self.getId()[2])  
391             else:             else:
392               out=data.interpolate(self.getFunctionSpace()).convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])               dat=data.interpolate(self.getFunctionSpace())
393               #out=data.interpolate(self.getFunctionSpace()).convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1],self.getId()[2])             id=self.getId()
394             if data.getRank()==0:             r=data.getRank()
395                return out[0]             if isinstance(id,list):
396                   out=[]
397                   for i in id:
398                      o=data.getValueOfGlobalDataPoint(*i)
399                      if data.getRank()==0:
400                         out.append(o[0])
401                      else:
402                         out.append(o)
403                   return out
404             else:             else:
405                return out               out=data.getValueOfGlobalDataPoint(*id)
406                 if data.getRank()==0:
407                    return out[0]
408                 else:
409                    return out
410          else:          else:
411             return data             return data
412    
413  # vim: expandtab shiftwidth=4:  class SolverSchemeException(Exception):
414       """
415       exceptions thrown by solvers
416       """
417       pass
418    
419    class IndefinitePreconditioner(SolverSchemeException):
420       """
421       the preconditioner is not positive definite.
422       """
423       pass
424    class MaxIterReached(SolverSchemeException):
425       """
426       maxium number of iteration steps is reached.
427       """
428       pass
429    class CorrectionFailed(SolverSchemeException):
430       """
431       no convergence has been achieved in the solution correction scheme.
432       """
433       pass
434    class IterationBreakDown(SolverSchemeException):
435       """
436       iteration scheme econouters an incurable breakdown.
437       """
438       pass
439    class NegativeNorm(SolverSchemeException):
440       """
441       a norm calculation returns a negative norm.
442       """
443       pass
444    
445    def PCG(b, Aprod, Msolve, bilinearform, atol=0, rtol=1.e-8, x=None, iter_max=100, initial_guess=True, verbose=False):
446       """
447       Solver for
448    
449       M{Ax=b}
450    
451       with a symmetric and positive definite operator A (more details required!).
452       It uses the conjugate gradient method with preconditioner M providing an approximation of A.
453    
454       The iteration is terminated if the C{stoppingcriterium} function return C{True}.
455    
456       For details on the preconditioned conjugate gradient method see the book:
457    
458       Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
459       T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
460       C. Romine, and H. van der Vorst.
461    
462       @param b: the right hand side of the liner system. C{b} is altered.
463       @type b: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
464       @param Aprod: returns the value Ax
465       @type Aprod: function C{Aprod(x)} where C{x} is of the same object like argument C{x}.
466           The returned object needs to be of the same type like argument C{b}.
467       @param Msolve: solves Mx=r
468       @type Msolve: function C{Msolve(r)} where C{r} is of the same type like argument C{b}.
469          The returned object needs to be of the same type like argument C{x}.
470       @param bilinearform: inner product C{<x,r>}
471       @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same type like argument C{x} and C{r} is.
472           The returned value is a C{float}.
473       @param stoppingcriterium: function which returns True if a stopping criterium is meet.
474           C{stoppingcriterium} has the arguments C{norm_r}, C{r} and C{x} giving the current
475           norm of the residual (=C{sqrt(bilinearform(Msolve(r),r)}), the current residual and
476           the current solution approximation. C{stoppingcriterium} is called in each iteration step.
477       @type stoppingcriterium: function that returns C{True} or C{False}
478       @param x: an initial guess for the solution. If no C{x} is given 0*b is used.
479       @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
480       @param iter_max: maximum number of iteration steps.
481       @type iter_max: C{int}
482       @return: the solution approximation and the corresponding residual
483       @rtype: C{tuple}
484       @warning: C{b} and C{x} are altered.
485       """
486       iter=0
487       if x==None:
488          x=0*b
489          r=b*1.
490          initial_guess=False
491       if initial_guess:
492          r=b+(-1)*Aprod(x)
493          x=x*1.
494       else:
495          x=x*0
496          r=b*1
497       rhat=Msolve(r)
498       d = rhat
499       rhat_dot_r = bilinearform(rhat, r)
500       if rhat_dot_r<0: raise NegativeNorm,"negative norm."
501       norm_b=math.sqrt(rhat_dot_r)
502       if verbose: print "PCG: initial residual norm = %e (absolute tolerance = %e)"%(norm_b, atol+rtol*norm_b)
503    
504       while not math.sqrt(rhat_dot_r) <= atol+rtol*norm_b:
505           iter+=1
506           if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
507    
508           q=Aprod(d)
509           alpha = rhat_dot_r / bilinearform(d, q)
510           x += alpha * d
511           r += (-alpha) * q
512    
513           rhat=Msolve(r)
514           rhat_dot_r_new = bilinearform(rhat, r)
515           beta = rhat_dot_r_new / rhat_dot_r
516           rhat+=beta * d
517           d=rhat
518    
519           rhat_dot_r = rhat_dot_r_new
520           if rhat_dot_r<0: raise NegativeNorm,"negative norm."
521           if verbose: print "PCG: iteration step %s: residual norm = %e"%(iter, math.sqrt(rhat_dot_r))
522       if verbose: print "PCG: tolerance reached after %s steps."%iter
523       return x,r
524    
525    class Defect(object):
526        """
527        defines a non-linear defect F(x) of a variable x
528        """
529        def __init__(self):
530            """
531            initialize defect
532            """
533            self.setDerivativeIncrementLength()
534    
535        def bilinearform(self, x0, x1):
536            """
537            returns the inner product of x0 and x1
538            @param x0: a value for x
539            @param x1: a value for x
540            @return: the inner product of x0 and x1
541            @rtype: C{float}
542            """
543            return 0
544          
545        def norm(self,x):
546            """
547            the norm of argument C{x}
548    
549            @param x: a value for x
550            @return: norm of argument x
551            @rtype: C{float}
552            @note: by default C{sqrt(self.bilinearform(x,x)} is retrurned.
553            """
554            s=self.bilinearform(x,x)
555            if s<0: raise NegativeNorm,"negative norm."
556            return math.sqrt(s)
557    
558    
559        def eval(self,x):
560            """
561            returns the value F of a given x
562    
563            @param x: value for which the defect C{F} is evalulated.
564            @return: value of the defect at C{x}
565            """
566            return 0
567    
568        def __call__(self,x):
569            return self.eval(x)
570    
571        def setDerivativeIncrementLength(self,inc=math.sqrt(util.EPSILON)):
572            """
573            sets the relative length of the increment used to approximate the derivative of the defect
574            the increment is inc*norm(x)/norm(v)*v in the direction of v with x as a staring point.
575    
576            @param inc: relative increment length
577            @type inc: positive C{float}
578            """
579            if inc<=0: raise ValueError,"positive increment required."
580            self.__inc=inc
581    
582        def getDerivativeIncrementLength(self):
583            """
584            returns the relative increment length used to approximate the derivative of the defect
585            @return: value of the defect at C{x}
586            @rtype: positive C{float}
587            """
588            return self.__inc
589    
590        def derivative(self, F0, x0, v, v_is_normalised=True):
591            """
592            returns the directional derivative at x0 in the direction of v
593    
594            @param F0: value of this defect at x0
595            @param x0: value at which derivative is calculated.
596            @param v: direction
597            @param v_is_normalised: is true to indicate that C{v} is nomalized (self.norm(v)=0)
598            @return: derivative of this defect at x0 in the direction of C{v}
599            @note: by default numerical evaluation (self.eval(x0+eps*v)-F0)/eps is used but this method
600            maybe oepsnew verwritten to use exact evalution.
601            """
602            normx=self.norm(x0)
603            if normx>0:
604                 epsnew = self.getDerivativeIncrementLength() * normx
605            else:
606                 epsnew = self.getDerivativeIncrementLength()
607            if not v_is_normalised:
608                normv=self.norm(v)
609                if normv<=0:
610                   return F0*0
611                else:
612                   epsnew /= normv
613            F1=self.eval(x0 + epsnew * v)
614            return (F1-F0)/epsnew
615    
616    ######################################    
617    def NewtonGMRES(defect, x, iter_max=100, sub_iter_max=20, atol=0,rtol=1.e-4, sub_tol_max=0.5, gamma=0.9, verbose=False):
618       """
619       solves a non-linear problem M{F(x)=0} for unknown M{x} using the stopping criterion:
620    
621       M{norm(F(x) <= atol + rtol * norm(F(x0)}
622      
623       where M{x0} is the initial guess.
624    
625       @param defect: object defining the the function M{F}, C{defect.norm} defines the M{norm} used in the stopping criterion.
626       @type defect: L{Defect}
627       @param x: initial guess for the solution, C{x} is altered.
628       @type x: any object type allowing basic operations such as  L{numarray.NumArray}, L{Data}
629       @param iter_max: maximum number of iteration steps
630       @type iter_max: positive C{int}
631       @param sub_iter_max:
632       @type sub_iter_max:
633       @param atol: absolute tolerance for the solution
634       @type atol: positive C{float}
635       @param rtol: relative tolerance for the solution
636       @type rtol: positive C{float}
637       @param gamma: tolerance safety factor for inner iteration
638       @type gamma: positive C{float}, less than 1
639       @param sub_tol_max: upper bound for inner tolerance.
640       @type sub_tol_max: positive C{float}, less than 1
641       @return: an approximation of the solution with the desired accuracy
642       @rtype: same type as the initial guess.
643       """
644       lmaxit=iter_max
645       if atol<0: raise ValueError,"atol needs to be non-negative."
646       if rtol<0: raise ValueError,"rtol needs to be non-negative."
647       if rtol+atol<=0: raise ValueError,"rtol or atol needs to be non-negative."
648       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
649       if sub_tol_max<=0 or sub_tol_max>=1: raise ValueError,"upper bound for inner tolerance for inner iteration (sub_tol_max =%s) needs to be positive and less than 1."%sub_tol_max
650    
651       F=defect(x)
652       fnrm=defect.norm(F)
653       stop_tol=atol + rtol*fnrm
654       sub_tol=sub_tol_max
655       if verbose: print "NewtonGMRES: initial residual = %e."%fnrm
656       if verbose: print "             tolerance = %e."%sub_tol
657       iter=1
658       #
659       # main iteration loop
660       #
661       while not fnrm<=stop_tol:
662                if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
663                #
664            #   adjust sub_tol_
665            #
666                if iter > 1:
667               rat=fnrm/fnrmo
668                   sub_tol_old=sub_tol
669               sub_tol=gamma*rat**2
670               if gamma*sub_tol_old**2 > .1: sub_tol=max(sub_tol,gamma*sub_tol_old**2)
671               sub_tol=max(min(sub_tol,sub_tol_max), .5*stop_tol/fnrm)
672            #
673            # calculate newton increment xc
674                #     if iter_max in __FDGMRES is reached MaxIterReached is thrown
675                #     if iter_restart -1 is returned as sub_iter
676                #     if  atol is reached sub_iter returns the numer of steps performed to get there
677                #
678                #  
679                if verbose: print "             subiteration (GMRES) is called with relative tolerance %e."%sub_tol
680                try:
681                   xc, sub_iter=__FDGMRES(F, defect, x, sub_tol*fnrm, iter_max=iter_max-iter, iter_restart=sub_iter_max)
682                except MaxIterReached:
683                   raise MaxIterReached,"maximum number of %s steps reached."%iter_max
684                if sub_iter<0:
685                   iter+=sub_iter_max
686                else:
687                   iter+=sub_iter
688                # ====
689            x+=xc
690                F=defect(x)
691            iter+=1
692                fnrmo, fnrm=fnrm, defect.norm(F)
693                if verbose: print "             step %s: residual %e."%(iter,fnrm)
694       if verbose: print "NewtonGMRES: completed after %s steps."%iter
695       return x
696    
697    def __givapp(c,s,vin):
698        """
699        apply a sequence of Givens rotations (c,s) to the recuirsively to the vector vin
700        @warning: C{vin} is altered.
701        """
702        vrot=vin
703        if isinstance(c,float):
704            vrot=[c*vrot[0]-s*vrot[1],s*vrot[0]+c*vrot[1]]
705        else:
706            for i in range(len(c)):
707                w1=c[i]*vrot[i]-s[i]*vrot[i+1]
708            w2=s[i]*vrot[i]+c[i]*vrot[i+1]
709                vrot[i:i+2]=w1,w2
710        return vrot
711    
712    def __FDGMRES(F0, defect, x0, atol, iter_max=100, iter_restart=20):
713       h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)
714       c=numarray.zeros(iter_restart,numarray.Float64)
715       s=numarray.zeros(iter_restart,numarray.Float64)
716       g=numarray.zeros(iter_restart,numarray.Float64)
717       v=[]
718    
719       rho=defect.norm(F0)
720       if rho<=0.: return x0*0
721      
722       v.append(-F0/rho)
723       g[0]=rho
724       iter=0
725       while rho > atol and iter<iter_restart-1:
726    
727        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
728    
729            p=defect.derivative(F0,x0,v[iter], v_is_normalised=True)
730        v.append(p)
731    
732        v_norm1=defect.norm(v[iter+1])
733    
734            # Modified Gram-Schmidt
735        for j in range(iter+1):
736             h[j,iter]=defect.bilinearform(v[j],v[iter+1])  
737             v[iter+1]-=h[j,iter]*v[j]
738          
739        h[iter+1,iter]=defect.norm(v[iter+1])
740        v_norm2=h[iter+1,iter]
741    
742            # Reorthogonalize if needed
743        if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
744            for j in range(iter+1):  
745               hr=defect.bilinearform(v[j],v[iter+1])
746                   h[j,iter]=h[j,iter]+hr
747                   v[iter+1] -= hr*v[j]
748    
749            v_norm2=defect.norm(v[iter+1])
750            h[iter+1,iter]=v_norm2
751            #   watch out for happy breakdown
752            if not v_norm2 == 0:
753                    v[iter+1]=v[iter+1]/h[iter+1,iter]
754    
755            #   Form and store the information for the new Givens rotation
756        if iter > 0 :
757            hhat=numarray.zeros(iter+1,numarray.Float64)
758            for i in range(iter+1) : hhat[i]=h[i,iter]
759            hhat=__givapp(c[0:iter],s[0:iter],hhat);
760                for i in range(iter+1) : h[i,iter]=hhat[i]
761    
762        mu=math.sqrt(h[iter,iter]*h[iter,iter]+h[iter+1,iter]*h[iter+1,iter])
763    
764        if mu!=0 :
765            c[iter]=h[iter,iter]/mu
766            s[iter]=-h[iter+1,iter]/mu
767            h[iter,iter]=c[iter]*h[iter,iter]-s[iter]*h[iter+1,iter]
768            h[iter+1,iter]=0.0
769            g[iter:iter+2]=__givapp(c[iter],s[iter],g[iter:iter+2])
770    
771            # Update the residual norm
772            rho=abs(g[iter+1])
773        iter+=1
774    
775       # At this point either iter > iter_max or rho < tol.
776       # It's time to compute x and leave.        
777       if iter > 0 :
778         y=numarray.zeros(iter,numarray.Float64)    
779         y[iter-1] = g[iter-1] / h[iter-1,iter-1]
780         if iter > 1 :  
781            i=iter-2  
782            while i>=0 :
783              y[i] = ( g[i] - numarray.dot(h[i,i+1:iter], y[i+1:iter])) / h[i,i]
784              i=i-1
785         xhat=v[iter-1]*y[iter-1]
786         for i in range(iter-1):
787        xhat += v[i]*y[i]
788       else :
789          xhat=v[0] * 0
790    
791       if iter<iter_restart-1:
792          stopped=iter
793       else:
794          stopped=-1
795    
796       return xhat,stopped
797    
798    ##############################################
799    def GMRES(b, Aprod, bilinearform, atol=0, rtol=1.e-8, x=None, iter_max=100, iter_restart=20, verbose=False):
800    ################################################
801       m=iter_restart
802       iter=0
803       if x == None:
804          xc=b*0
805       else:
806          xc=x*1
807       if rtol>0:
808          r_dot_r = bilinearform(b, b)
809          if r_dot_r<0: raise NegativeNorm,"negative norm."
810          atol2=atol+rtol*math.sqrt(r_dot_r)
811          if verbose: print "GMRES: norm of right hand side = %e (absolute tolerance = %e)"%(math.sqrt(r_dot_r), atol2)
812       else:
813          atol2=atol
814          if verbose: print "GMRES: absolute tolerance = %e"%atol2
815      
816       while True:
817          if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached"%iter_max
818          xc,stopped=__GMRESm(b, Aprod, bilinearform, atol2, x=xc, iter_max=iter_max-iter, iter_restart=m, verbose=verbose)
819          iter+=iter_restart    
820          if stopped: break
821          if verbose: print "GMRES: restart."
822       if verbose: print "GMRES: tolerance has reached."
823       return xc
824    
825    def __GMRESm(b, Aprod, bilinearform, atol, x, iter_max=100, iter_restart=20, verbose=False):
826       iter=0
827      
828       h=numarray.zeros((iter_restart+1,iter_restart),numarray.Float64)
829       c=numarray.zeros(iter_restart,numarray.Float64)
830       s=numarray.zeros(iter_restart,numarray.Float64)
831       g=numarray.zeros(iter_restart+1,numarray.Float64)
832       v=[]
833    
834       r=b-Aprod(x)
835       r_dot_r = bilinearform(r, r)
836       if r_dot_r<0: raise NegativeNorm,"negative norm."
837       rho=math.sqrt(r_dot_r)
838      
839       v.append(r/rho)
840       g[0]=rho
841    
842       if verbose: print "GMRES: initial residual %e (absolute tolerance = %e)"%(rho,atol)
843       while not (rho<=atol or iter==iter_restart):
844    
845        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
846    
847        p=Aprod(v[iter])
848        v.append(p)
849    
850        v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
851    
852    # Modified Gram-Schmidt
853        for j in range(iter+1):
854          h[j,iter]=bilinearform(v[j],v[iter+1])  
855          v[iter+1]-=h[j,iter]*v[j]
856          
857        h[iter+1,iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
858        v_norm2=h[iter+1,iter]
859    
860    # Reorthogonalize if needed
861        if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
862         for j in range(iter+1):  
863            hr=bilinearform(v[j],v[iter+1])
864                h[j,iter]=h[j,iter]+hr
865                v[iter+1] -= hr*v[j]
866    
867         v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
868         h[iter+1,iter]=v_norm2
869    
870    #   watch out for happy breakdown
871            if not v_norm2 == 0:
872             v[iter+1]=v[iter+1]/h[iter+1,iter]
873    
874    #   Form and store the information for the new Givens rotation
875        if iter > 0: h[:iter+1,iter]=__givapp(c[:iter],s[:iter],h[:iter+1,iter])
876        mu=math.sqrt(h[iter,iter]*h[iter,iter]+h[iter+1,iter]*h[iter+1,iter])
877    
878        if mu!=0 :
879            c[iter]=h[iter,iter]/mu
880            s[iter]=-h[iter+1,iter]/mu
881            h[iter,iter]=c[iter]*h[iter,iter]-s[iter]*h[iter+1,iter]
882            h[iter+1,iter]=0.0
883            g[iter:iter+2]=__givapp(c[iter],s[iter],g[iter:iter+2])
884    # Update the residual norm
885                  
886            rho=abs(g[iter+1])
887            if verbose: print "GMRES: iteration step %s: residual %e"%(iter,rho)
888        iter+=1
889    
890    # At this point either iter > iter_max or rho < tol.
891    # It's time to compute x and leave.        
892    
893       if verbose: print "GMRES: iteration stopped after %s step."%iter
894       if iter > 0 :
895         y=numarray.zeros(iter,numarray.Float64)    
896         y[iter-1] = g[iter-1] / h[iter-1,iter-1]
897         if iter > 1 :  
898            i=iter-2  
899            while i>=0 :
900              y[i] = ( g[i] - numarray.dot(h[i,i+1:iter], y[i+1:iter])) / h[i,i]
901              i=i-1
902         xhat=v[iter-1]*y[iter-1]
903         for i in range(iter-1):
904        xhat += v[i]*y[i]
905       else:
906         xhat=v[0] * 0
907    
908       x += xhat
909       if iter<iter_restart-1:
910          stopped=True
911       else:
912          stopped=False
913    
914       return x,stopped
915    
916    #################################################
917    def MINRES(b, Aprod, Msolve, bilinearform, atol=0, rtol=1.e-8, x=None, iter_max=100):
918    #################################################
919        #
920        #  minres solves the system of linear equations Ax = b
921        #  where A is a symmetric matrix (possibly indefinite or singular)
922        #  and b is a given vector.
923        #  
924        #  "A" may be a dense or sparse matrix (preferably sparse!)
925        #  or the name of a function such that
926        #               y = A(x)
927        #  returns the product y = Ax for any given vector x.
928        #
929        #  "M" defines a positive-definite preconditioner M = C C'.
930        #  "M" may be a dense or sparse matrix (preferably sparse!)
931        #  or the name of a function such that
932        #  solves the system My = x for any given vector x.
933        #
934        #
935        
936        #------------------------------------------------------------------
937        # Set up y and v for the first Lanczos vector v1.
938        # y  =  beta1 P' v1,  where  P = C**(-1).
939        # v is really P' v1.
940        #------------------------------------------------------------------
941        if x==None:
942          x=0*b
943        else:
944          b += (-1)*Aprod(x)
945    
946        r1    = b
947        y = Msolve(b)
948        beta1 = bilinearform(y,b)
949    
950        if beta1< 0: raise NegativeNorm,"negative norm."
951    
952        #  If b = 0 exactly, stop with x = 0.
953        if beta1==0: return x*0.
954    
955        if beta1> 0:
956          beta1  = math.sqrt(beta1)      
957    
958        #------------------------------------------------------------------
959        # Initialize quantities.
960        # ------------------------------------------------------------------
961        iter   = 0
962        Anorm = 0
963        ynorm = 0
964        oldb   = 0
965        beta   = beta1
966        dbar   = 0
967        epsln  = 0
968        phibar = beta1
969        rhs1   = beta1
970        rhs2   = 0
971        rnorm  = phibar
972        tnorm2 = 0
973        ynorm2 = 0
974        cs     = -1
975        sn     = 0
976        w      = b*0.
977        w2     = b*0.
978        r2     = r1
979        eps    = 0.0001
980    
981        #---------------------------------------------------------------------
982        # Main iteration loop.
983        # --------------------------------------------------------------------
984        while not rnorm<=atol+rtol*Anorm*ynorm:    #  checks ||r|| < (||A|| ||x||) * TOL
985    
986        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
987            iter    = iter  +  1
988    
989            #-----------------------------------------------------------------
990            # Obtain quantities for the next Lanczos vector vk+1, k = 1, 2,...
991            # The general iteration is similar to the case k = 1 with v0 = 0:
992            #
993            #   p1      = Operator * v1  -  beta1 * v0,
994            #   alpha1  = v1'p1,
995            #   q2      = p2  -  alpha1 * v1,
996            #   beta2^2 = q2'q2,
997            #   v2      = (1/beta2) q2.
998            #
999            # Again, y = betak P vk,  where  P = C**(-1).
1000            #-----------------------------------------------------------------
1001            s = 1/beta                 # Normalize previous vector (in y).
1002            v = s*y                    # v = vk if P = I
1003        
1004            y      = Aprod(v)
1005        
1006            if iter >= 2:
1007              y = y - (beta/oldb)*r1
1008    
1009            alfa   = bilinearform(v,y)              # alphak
1010            y      += (- alfa/beta)*r2
1011            r1     = r2
1012            r2     = y
1013            y = Msolve(r2)
1014            oldb   = beta                           # oldb = betak
1015            beta   = bilinearform(y,r2)             # beta = betak+1^2
1016            if beta < 0: raise NegativeNorm,"negative norm."
1017    
1018            beta   = math.sqrt( beta )
1019            tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta
1020            
1021            if iter==1:                 # Initialize a few things.
1022              gmax   = abs( alfa )      # alpha1
1023              gmin   = gmax             # alpha1
1024    
1025            # Apply previous rotation Qk-1 to get
1026            #   [deltak epslnk+1] = [cs  sn][dbark    0   ]
1027            #   [gbar k dbar k+1]   [sn -cs][alfak betak+1].
1028        
1029            oldeps = epsln
1030            delta  = cs * dbar  +  sn * alfa  # delta1 = 0         deltak
1031            gbar   = sn * dbar  -  cs * alfa  # gbar 1 = alfa1     gbar k
1032            epsln  =               sn * beta  # epsln2 = 0         epslnk+1
1033            dbar   =            -  cs * beta  # dbar 2 = beta2     dbar k+1
1034    
1035            # Compute the next plane rotation Qk
1036    
1037            gamma  = math.sqrt(gbar*gbar+beta*beta)  # gammak
1038            gamma  = max(gamma,eps)
1039            cs     = gbar / gamma             # ck
1040            sn     = beta / gamma             # sk
1041            phi    = cs * phibar              # phik
1042            phibar = sn * phibar              # phibark+1
1043    
1044            # Update  x.
1045    
1046            denom = 1/gamma
1047            w1    = w2
1048            w2    = w
1049            w     = (v - oldeps*w1 - delta*w2) * denom
1050            x     +=  phi*w
1051    
1052            # Go round again.
1053    
1054            gmax   = max(gmax,gamma)
1055            gmin   = min(gmin,gamma)
1056            z      = rhs1 / gamma
1057            ynorm2 = z*z  +  ynorm2
1058            rhs1   = rhs2 -  delta*z
1059            rhs2   =      -  epsln*z
1060    
1061            # Estimate various norms and test for convergence.
1062    
1063            Anorm  = math.sqrt( tnorm2 )
1064            ynorm  = math.sqrt( ynorm2 )
1065    
1066            rnorm  = phibar
1067    
1068        return x
1069    
1070    def TFQMR(b, Aprod, Msolve, bilinearform, atol=0, rtol=1.e-8, x=None, iter_max=100):
1071    
1072    # TFQMR solver for linear systems
1073    #
1074    #
1075    # initialization
1076    #
1077      errtol = math.sqrt(bilinearform(b,b))
1078      norm_b=errtol
1079      kmax  = iter_max
1080      error = []
1081    
1082      if math.sqrt(bilinearform(x,x)) != 0.0:
1083        r = b - Aprod(x)
1084      else:
1085        r = b
1086    
1087      r=Msolve(r)
1088    
1089      u1=0
1090      u2=0
1091      y1=0
1092      y2=0
1093    
1094      w = r
1095      y1 = r
1096      iter = 0
1097      d = 0
1098      
1099      v = Msolve(Aprod(y1))
1100      u1 = v
1101      
1102      theta = 0.0;
1103      eta = 0.0;
1104      tau = math.sqrt(bilinearform(r,r))
1105      error = [ error, tau ]
1106      rho = tau * tau
1107    #
1108    #  TFQMR iteration
1109    #
1110    #  while ( iter < kmax-1 ):
1111      
1112      while not tau<=atol+rtol*norm_b:
1113        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
1114    
1115        sigma = bilinearform(r,v)
1116    
1117        if ( sigma == 0.0 ):
1118          raise 'TFQMR breakdown, sigma=0'
1119        
1120    
1121        alpha = rho / sigma
1122    
1123        for j in range(2):
1124    #
1125    #   Compute y2 and u2 only if you have to
1126    #
1127          if ( j == 1 ):
1128            y2 = y1 - alpha * v
1129            u2 = Msolve(Aprod(y2))
1130    
1131          m = 2 * (iter+1) - 2 + (j+1)
1132          if j==0:
1133             w = w - alpha * u1
1134             d = y1 + ( theta * theta * eta / alpha ) * d
1135          if j==1:
1136             w = w - alpha * u2
1137             d = y2 + ( theta * theta * eta / alpha ) * d
1138    
1139          theta = math.sqrt(bilinearform(w,w))/ tau
1140          c = 1.0 / math.sqrt ( 1.0 + theta * theta )
1141          tau = tau * theta * c
1142          eta = c * c * alpha
1143          x = x + eta * d
1144    #
1145    #  Try to terminate the iteration at each pass through the loop
1146    #
1147         # if ( tau * math.sqrt ( m + 1 ) <= errtol ):
1148         #   error = [ error, tau ]
1149         #   total_iters = iter
1150         #   break
1151          
1152    
1153        if ( rho == 0.0 ):
1154          raise 'TFQMR breakdown, rho=0'
1155        
1156    
1157        rhon = bilinearform(r,w)
1158        beta = rhon / rho;
1159        rho = rhon;
1160        y1 = w + beta * y2;
1161        u1 = Msolve(Aprod(y1))
1162        v = u1 + beta * ( u2 + beta * v )
1163        error = [ error, tau ]
1164        total_iters = iter
1165        
1166        iter = iter + 1
1167    
1168      return x
1169    
1170    
1171    #############################################
1172    
1173    class ArithmeticTuple(object):
1174       """
1175       tuple supporting inplace update x+=y and scaling x=a*y where x,y is an ArithmeticTuple and a is a float.
1176    
1177       example of usage:
1178    
1179       from esys.escript import Data
1180       from numarray import array
1181       a=Data(...)
1182       b=array([1.,4.])
1183       x=ArithmeticTuple(a,b)
1184       y=5.*x
1185    
1186       """
1187       def __init__(self,*args):
1188           """
1189           initialize object with elements args.
1190    
1191           @param args: tuple of object that support implace add (x+=y) and scaling (x=a*y)
1192           """
1193           self.__items=list(args)
1194    
1195       def __len__(self):
1196           """
1197           number of items
1198    
1199           @return: number of items
1200           @rtype: C{int}
1201           """
1202           return len(self.__items)
1203    
1204       def __getitem__(self,index):
1205           """
1206           get an item
1207    
1208           @param index: item to be returned
1209           @type index: C{int}
1210           @return: item with index C{index}
1211           """
1212           return self.__items.__getitem__(index)
1213    
1214       def __mul__(self,other):
1215           """
1216           scaling from the right
1217    
1218           @param other: scaling factor
1219           @type other: C{float}
1220           @return: itemwise self*other
1221           @rtype: L{ArithmeticTuple}
1222           """
1223           out=[]
1224           try:  
1225               l=len(other)
1226               if l!=len(self):
1227                   raise ValueError,"length of of arguments don't match."
1228               for i in range(l): out.append(self[i]*other[i])
1229           except TypeError:
1230           for i in range(len(self)): out.append(self[i]*other)
1231           return ArithmeticTuple(*tuple(out))
1232    
1233       def __rmul__(self,other):
1234           """
1235           scaling from the left
1236    
1237           @param other: scaling factor
1238           @type other: C{float}
1239           @return: itemwise other*self
1240           @rtype: L{ArithmeticTuple}
1241           """
1242           out=[]
1243           try:  
1244               l=len(other)
1245               if l!=len(self):
1246                   raise ValueError,"length of of arguments don't match."
1247               for i in range(l): out.append(other[i]*self[i])
1248           except TypeError:
1249           for i in range(len(self)): out.append(other*self[i])
1250           return ArithmeticTuple(*tuple(out))
1251    
1252       def __div__(self,other):
1253           """
1254           dividing from the right
1255    
1256           @param other: scaling factor
1257           @type other: C{float}
1258           @return: itemwise self/other
1259           @rtype: L{ArithmeticTuple}
1260           """
1261           return self*(1/other)
1262    
1263       def __rdiv__(self,other):
1264           """
1265           dividing from the left
1266    
1267           @param other: scaling factor
1268           @type other: C{float}
1269           @return: itemwise other/self
1270           @rtype: L{ArithmeticTuple}
1271           """
1272           out=[]
1273           try:  
1274               l=len(other)
1275               if l!=len(self):
1276                   raise ValueError,"length of of arguments don't match."
1277               for i in range(l): out.append(other[i]/self[i])
1278           except TypeError:
1279           for i in range(len(self)): out.append(other/self[i])
1280           return ArithmeticTuple(*tuple(out))
1281      
1282       def __iadd__(self,other):
1283           """
1284           in-place add of other to self
1285    
1286           @param other: increment
1287           @type other: C{ArithmeticTuple}
1288           """
1289           if len(self) != len(other):
1290               raise ValueError,"tuple length must match."
1291           for i in range(len(self)):
1292               self.__items[i]+=other[i]
1293           return self
1294    
1295       def __add__(self,other):
1296           """
1297           add other to self
1298    
1299           @param other: increment
1300           @type other: C{ArithmeticTuple}
1301           """
1302           out=[]
1303           try:  
1304               l=len(other)
1305               if l!=len(self):
1306                   raise ValueError,"length of of arguments don't match."
1307               for i in range(l): out.append(self[i]+other[i])
1308           except TypeError:
1309           for i in range(len(self)): out.append(self[i]+other)
1310           return ArithmeticTuple(*tuple(out))
1311    
1312       def __sub__(self,other):
1313           """
1314           subtract other from self
1315    
1316           @param other: increment
1317           @type other: C{ArithmeticTuple}
1318           """
1319           out=[]
1320           try:  
1321               l=len(other)
1322               if l!=len(self):
1323                   raise ValueError,"length of of arguments don't match."
1324               for i in range(l): out.append(self[i]-other[i])
1325           except TypeError:
1326           for i in range(len(self)): out.append(self[i]-other)
1327           return ArithmeticTuple(*tuple(out))
1328      
1329       def __isub__(self,other):
1330           """
1331           subtract other from self
1332    
1333           @param other: increment
1334           @type other: C{ArithmeticTuple}
1335           """
1336           if len(self) != len(other):
1337               raise ValueError,"tuple length must match."
1338           for i in range(len(self)):
1339               self.__items[i]-=other[i]
1340           return self
1341    
1342       def __neg__(self):
1343           """
1344           negate
1345    
1346           """
1347           out=[]
1348           for i in range(len(self)):
1349               out.append(-self[i])
1350           return ArithmeticTuple(*tuple(out))
1351    
1352    
1353    class HomogeneousSaddlePointProblem(object):
1354          """
1355          This provides a framwork for solving linear homogeneous saddle point problem of the form
1356    
1357                 Av+B^*p=f
1358                 Bv    =0
1359    
1360          for the unknowns v and p and given operators A and B and given right hand side f.
1361          B^* is the adjoint operator of B.
1362          """
1363          def __init__(self,**kwargs):
1364            self.setTolerance()
1365            self.setAbsoluteTolerance()
1366            self.setSubToleranceReductionFactor()
1367    
1368          #=============================================================
1369          def initialize(self):
1370            """
1371            initialize the problem (overwrite)
1372            """
1373            pass
1374    
1375          def B(self,v):
1376            """
1377            returns Bv (overwrite)
1378            @rtype: equal to the type of p
1379    
1380            @note: boundary conditions on p should be zero!
1381            """
1382            raise NotImplementedError,"no operator B implemented."
1383    
1384          def inner_pBv(self,p,Bv):
1385             """
1386             returns inner product of element p and Bv  (overwrite)
1387            
1388             @type p: equal to the type of p
1389             @type Bv: equal to the type of result of operator B
1390             @rtype: C{float}
1391    
1392             @rtype: equal to the type of p
1393             """
1394             raise NotImplementedError,"no inner product for p implemented."
1395    
1396          def inner_p(self,p0,p1):
1397             """
1398             returns inner product of element p0 and p1  (overwrite)
1399            
1400             @type p0: equal to the type of p
1401             @type p1: equal to the type of p
1402             @rtype: equal to the type of p
1403             """
1404             raise NotImplementedError,"no inner product for p implemented."
1405    
1406          def inner_v(self,v0,v1):
1407             """
1408             returns inner product of two element v0 and v1  (overwrite)
1409            
1410             @type v0: equal to the type of v
1411             @type v1: equal to the type of v
1412             @rtype: C{float}
1413    
1414             @rtype: equal to the type of v
1415             """
1416             raise NotImplementedError,"no inner product for v implemented."
1417             pass
1418    
1419          def solve_A(self,u,p):
1420             """
1421             solves Av=f-Au-B^*p with accuracy self.getSubProblemTolerance() (overwrite)
1422    
1423             @rtype: equal to the type of v
1424             @note: boundary conditions on v should be zero!
1425             """
1426             raise NotImplementedError,"no operator A implemented."
1427    
1428          def solve_prec(self,p):
1429             """
1430             provides a preconditioner for BA^{-1}B^* with accuracy self.getSubProblemTolerance() (overwrite)
1431    
1432             @rtype: equal to the type of p
1433             @note: boundary conditions on p should be zero!
1434             """
1435             raise NotImplementedError,"no preconditioner for Schur complement implemented."
1436          #=============================================================
1437          def __Aprod_PCG(self,p):
1438              # return (v,Bv) with v=A^-1B*p
1439              #solve Av =B^*p as Av =f-Az-B^*(-p)
1440              v=self.solve_A(self.__z,-p)
1441              return ArithmeticTuple(v, self.B(v))
1442    
1443          def __inner_PCG(self,p,r):
1444             return self.inner_pBv(p,r[1])
1445    
1446          def __Msolve_PCG(self,r):
1447              return self.solve_prec(r[1])
1448          #=============================================================
1449          def __Aprod_GMRES(self,x):
1450              # return w,q  from v, p
1451              # solve Aw =Av+B^*p as Aw =f-A(z-v)-B^*(-p)
1452              #  and  Sq=B(v-w)
1453              v=x[0]
1454              p=x[1]
1455              w=self.solve_A(self.__z-v,-p)
1456              Bw=self.B(w-v)
1457          q=self.solve_prec(Bw)
1458              return ArithmeticTuple(w,q)
1459    
1460          def __inner_GMRES(self,a1,a2):
1461             return self.inner_p(a1[1],a2[1])+self.inner_v(a1[0],a2[0])
1462    
1463          #=============================================================
1464          def norm(self,v,p):
1465            f=self.inner_p(p,p)+self.inner_v(v,v)
1466            if f<0:
1467                raise ValueError,"negative norm."
1468            return math.sqrt(f)
1469    
1470          def solve(self,v,p,max_iter=20, verbose=False, show_details=False, useUzawa=True, iter_restart=20,max_correction_steps=3):
1471             """
1472             solves the saddle point problem using initial guesses v and p.
1473    
1474             @param v: initial guess for velocity
1475             @param p: initial guess for pressure
1476             @type v: L{Data}
1477             @type p: L{Data}
1478             @param useUzawa: indicate the usage of the Uzawa scheme. Otherwise the problem is solved in coupled form.
1479             @param max_iter: maximum number of iteration steps per correction attempt.
1480             @param verbose: show information on the progress of the saddlepoint problem solver.
1481             @param show_details: show details of the sub problems solves
1482             @param iter_restart: restart the iteration after C{iter_restart} steps (only used if useUzaw=False)
1483             @param max_correction_steps: maximum number of iteration steps in the attempt get |Bv| to zero.
1484             @return: new approximations for velocity and pressure
1485             @type useUzawa: C{bool}
1486             @type max_iter: C{int}
1487             @type verbose: C{bool}
1488             @type show_details: C{bool}
1489             @type iter_restart: C{int}
1490             @type max_correction_steps: C{int}
1491             @rtype: C{tuple} of L{Data} objects
1492             """
1493             self.verbose=verbose
1494             self.show_details=show_details and self.verbose
1495             #=================================================================================
1496             # Az=f is solved as A(z-v)=f-Av-B^*0 (typically with z-v=0 on boundary)
1497             self.setSubProblemTolerance(self.getTolerance(), apply_reduction=True)
1498             self.__z=v+self.solve_A(v,p*0)
1499             Bz=self.B(self.__z)
1500             p0=self.solve_prec(Bz)
1501             # tolerances:
1502             rtol=self.getTolerance()
1503             atol=self.getAbsoluteTolerance()
1504             f0=self.norm(self.__z,p0)
1505             if not f0>0:
1506                return self.__z, p*0
1507             ATOL=rtol*f0+atol
1508             if not ATOL>0:
1509                raise ValueError,"overall absolute tolerance needs to be positive."
1510             if self.verbose: print "saddle point solver: initial residual %e, tolerance = %e."%(f0,ATOL)
1511             # initialization
1512             self.iter=0
1513             correction_step=0
1514             converged=False
1515             # initial guess:
1516             q=p*1
1517             if not useUzawa:
1518               u=v
1519             while not converged :
1520                if useUzawa:
1521                   self.setSubProblemTolerance(ATOL/f0, apply_reduction=True)
1522                   # assume p is known: then v=z-A^-1B^*p
1523                   #      
1524                   # which leads to BA^-1B^*p = Bz
1525                   #
1526                   # note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bu
1527                   #
1528                   # note that q=BA^-1B^*p = B (A^-1B^*p)=Bw  with Aw=B^*p = f - Az - B^*(-p)
1529                   #
1530                   if self.verbose: print "enter PCG method (iter_max=%s, atol=%e, subtolerance=%e)"%(max_iter,ATOL, self.getSubProblemTolerance())
1531                   q,r=PCG(ArithmeticTuple(self.__z,Bz),self.__Aprod_PCG,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL, rtol=0.,iter_max=max_iter, x=q, verbose=self.verbose)
1532               u=r[0]  
1533               Bu=r[1]  
1534                else:
1535                   self.setSubProblemTolerance(ATOL/f0, apply_reduction=True)
1536                   #
1537                   #  with v=dv+z
1538                   #
1539                   #   A dv + B p = A (v-z) + Bp = 0
1540                   #   B dv       =  - Bz
1541                   #
1542                   # apply the preconditioner [[A^{-1} 0][(S^{-1} B A^{-1})  -S^{-1}]] to the problem
1543                   # the right hand side is then [0, S^{-1} Bz ]
1544                   #
1545                   #          
1546                   if self.verbose: print "enter GMRES (iter_max=%s, atol=%e, subtolerance=%e)"%(max_iter,ATOL,self.getSubProblemTolerance())
1547                   x=GMRES(ArithmeticTuple(0*self.__z,p0),self.__Aprod_GMRES,self.__inner_GMRES,atol=ATOL, rtol=0.,iter_max=max_iter, x=ArithmeticTuple(u-self.__z,q),iter_restart=iter_restart, verbose=self.verbose)
1548                   self.setSubProblemTolerance(ATOL/f0, apply_reduction=False)
1549               u=self.__z+x[0]
1550                   q=x[1]
1551               Bu=self.B(u)
1552                # now we check if |Bu| ~ 0 or more precise |Bu|  <= rtol * |v|
1553                nu=self.inner_v(u,u)
1554                nz=self.inner_v(self.__z,self.__z)
1555                p2=self.solve_prec(Bu)
1556                nBu=self.inner_p(p2,p2)
1557                if not nu>=0 and not nBu>=0: raise NegativeNorm,"negative norm."
1558                nu=math.sqrt(nu)
1559                nBu=math.sqrt(nBu)
1560                if self.verbose: print "saddle point solver: norm v= %e (Bv = %e)"%(nu,nBu)
1561                QTOL=atol+nu*rtol
1562                if nBu <= QTOL:
1563                    converged=True
1564                else:
1565                    ATOL=QTOL/nBu*ATOL*0.3
1566                    if self.verbose: print "correction step %s: tolerance reduced to %e."%(correction_step,ATOL)
1567                    converged=False
1568                correction_step+=1
1569                if correction_step>max_correction_steps:
1570                   raise CorrectionFailed,"Given up after %d correction steps."%correction_step
1571             if self.verbose: print "saddle point solver: tolerance reached."
1572         return u,q
1573    
1574          #==========================================================================================================================
1575          def setTolerance(self,tolerance=1.e-4):
1576             """
1577             sets the relative tolerance for (v,p)
1578    
1579             @param tolerance: tolerance to be used
1580             @type tolerance: non-negative C{float}
1581             """
1582             if tolerance<0:
1583                 raise ValueError,"tolerance must be positive."
1584             self.__rtol=tolerance
1585          def getTolerance(self):
1586             """
1587             returns the relative tolerance
1588    
1589             @return: relative tolerance
1590             @rtype: C{float}
1591             """
1592             return self.__rtol
1593          def setAbsoluteTolerance(self,tolerance=0.):
1594             """
1595             sets the absolute tolerance
1596    
1597             @param tolerance: tolerance to be used
1598             @type tolerance: non-negative C{float}
1599             """
1600             if tolerance<0:
1601                 raise ValueError,"tolerance must be non-negative."
1602             self.__atol=tolerance
1603          def getAbsoluteTolerance(self):
1604             """
1605             returns the absolute tolerance
1606    
1607             @return: absolute tolerance
1608             @rtype: C{float}
1609             """
1610             return self.__atol
1611          def setSubToleranceReductionFactor(self,reduction=None):
1612             """
1613             sets the reduction factor for the tolerance used to solve the subproblem in each iteration step.
1614             if set to C{None} the square of the problem tolerance is used.
1615        
1616             @param reduction: reduction factor.
1617             @type reduction: positive C{float} less or equal one or C{None}
1618             """
1619             if not reduction == None:
1620                  if reduction<=0 or reduction>1:
1621                      raise ValueError,"reduction factor (=%e) must be positive and less or equal one."%reduction
1622             self.__reduction=reduction
1623          def getSubToleranceReductionFactor(self):
1624             """
1625             returns the subproblem reduction factor
1626    
1627             @return: subproblem reduction factor
1628             @rtype: C{float} or C{None}
1629             """
1630             return self.__reduction
1631    
1632          def setSubProblemTolerance(self,rtol,apply_reduction=False):
1633             """
1634             sets the relative tolerance to solve the subproblem(s).
1635    
1636             @param rtol: relative tolerence
1637             @type rtol: positive C{float}
1638             @param apply_reduction: if set the the reduction factor is applied.
1639             @type apply_reduction: C{bool}
1640             """
1641             if rtol<=0:
1642                 raise ValueError,"tolerance must be positive."
1643             if apply_reduction:
1644                 if self.__reduction == None:
1645                     self.__sub_tol=max(rtol**2,util.EPSILON*10)
1646                 else:
1647                     self.__sub_tol=max(self.__reduction*rtol,util.EPSILON*10)
1648             else:
1649                 self.__sub_tol=rtol
1650          def getSubProblemTolerance(self):
1651             """
1652             returns the subproblem reduction factor
1653    
1654             @return: subproblem reduction factor
1655             @rtype: C{float}
1656             """
1657             return self.__sub_tol
1658    
1659    def MaskFromBoundaryTag(domain,*tags):
1660       """
1661       create a mask on the Solution(domain) function space which one for samples
1662       that touch the boundary tagged by tags.
1663    
1664       usage: m=MaskFromBoundaryTag(domain,"left", "right")
1665    
1666       @param domain: a given domain
1667       @type domain: L{escript.Domain}
1668       @param tags: boundray tags
1669       @type tags: C{str}
1670       @return: a mask which marks samples that are touching the boundary tagged by any of the given tags.
1671       @rtype: L{escript.Data} of rank 0
1672       """
1673       pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1674       d=escript.Scalar(0.,escript.FunctionOnBoundary(domain))
1675       for t in tags: d.setTaggedValue(t,1.)
1676       pde.setValue(y=d)
1677       return util.whereNonZero(pde.getRightHandSide())
1678    #============================================================================================================================================
1679    class SaddlePointProblem(object):
1680       """
1681       This implements a solver for a saddlepoint problem
1682    
1683       M{f(u,p)=0}
1684       M{g(u)=0}
1685    
1686       for u and p. The problem is solved with an inexact Uszawa scheme for p:
1687    
1688       M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})}
1689       M{Q_g (p^{k+1}-p^{k}) =   g(u^{k+1})}
1690    
1691       where Q_f is an approximation of the Jacobiean A_f of f with respect to u  and Q_f is an approximation of
1692       A_g A_f^{-1} A_g with A_g is the jacobiean of g with respect to p. As a the construction of a 'proper'
1693       Q_g can be difficult, non-linear conjugate gradient method is applied to solve for p, so Q_g plays
1694       in fact the role of a preconditioner.
1695       """
1696       def __init__(self,verbose=False,*args):
1697           """
1698           initializes the problem
1699    
1700           @param verbose: switches on the printing out some information
1701           @type verbose: C{bool}
1702           @note: this method may be overwritten by a particular saddle point problem
1703           """
1704           print "SaddlePointProblem should not be used anymore!"
1705           if not isinstance(verbose,bool):
1706                raise TypeError("verbose needs to be of type bool.")
1707           self.__verbose=verbose
1708           self.relaxation=1.
1709           DeprecationWarning("SaddlePointProblem should not be used anymore.")
1710    
1711       def trace(self,text):
1712           """
1713           prints text if verbose has been set
1714    
1715           @param text: a text message
1716           @type text: C{str}
1717           """
1718           if self.__verbose: print "%s: %s"%(str(self),text)
1719    
1720       def solve_f(self,u,p,tol=1.e-8):
1721           """
1722           solves
1723    
1724           A_f du = f(u,p)
1725    
1726           with tolerance C{tol} and return du. A_f is Jacobiean of f with respect to u.
1727    
1728           @param u: current approximation of u
1729           @type u: L{escript.Data}
1730           @param p: current approximation of p
1731           @type p: L{escript.Data}
1732           @param tol: tolerance expected for du
1733           @type tol: C{float}
1734           @return: increment du
1735           @rtype: L{escript.Data}
1736           @note: this method has to be overwritten by a particular saddle point problem
1737           """
1738           pass
1739    
1740       def solve_g(self,u,tol=1.e-8):
1741           """
1742           solves
1743    
1744           Q_g dp = g(u)
1745    
1746           with Q_g is a preconditioner for A_g A_f^{-1} A_g with  A_g is the jacobiean of g with respect to p.
1747    
1748           @param u: current approximation of u
1749           @type u: L{escript.Data}
1750           @param tol: tolerance expected for dp
1751           @type tol: C{float}
1752           @return: increment dp
1753           @rtype: L{escript.Data}
1754           @note: this method has to be overwritten by a particular saddle point problem
1755           """
1756           pass
1757    
1758       def inner(self,p0,p1):
1759           """
1760           inner product of p0 and p1 approximating p. Typically this returns integrate(p0*p1)
1761           @return: inner product of p0 and p1
1762           @rtype: C{float}
1763           """
1764           pass
1765    
1766       subiter_max=3
1767       def solve(self,u0,p0,tolerance=1.e-6,tolerance_u=None,iter_max=100,accepted_reduction=0.995,relaxation=None):
1768            """
1769            runs the solver
1770    
1771            @param u0: initial guess for C{u}
1772            @type u0: L{esys.escript.Data}
1773            @param p0: initial guess for C{p}
1774            @type p0: L{esys.escript.Data}
1775            @param tolerance: tolerance for relative error in C{u} and C{p}
1776            @type tolerance: positive C{float}
1777            @param tolerance_u: tolerance for relative error in C{u} if different from C{tolerance}
1778            @type tolerance_u: positive C{float}
1779            @param iter_max: maximum number of iteration steps.
1780            @type iter_max: C{int}
1781            @param accepted_reduction: if the norm  g cannot be reduced by C{accepted_reduction} backtracking to adapt the
1782                                       relaxation factor. If C{accepted_reduction=None} no backtracking is used.
1783            @type accepted_reduction: positive C{float} or C{None}
1784            @param relaxation: initial relaxation factor. If C{relaxation==None}, the last relaxation factor is used.
1785            @type relaxation: C{float} or C{None}
1786            """
1787            tol=1.e-2
1788            if tolerance_u==None: tolerance_u=tolerance
1789            if not relaxation==None: self.relaxation=relaxation
1790            if accepted_reduction ==None:
1791                  angle_limit=0.
1792            elif accepted_reduction>=1.:
1793                  angle_limit=0.
1794            else:
1795                  angle_limit=util.sqrt(1-accepted_reduction**2)
1796            self.iter=0
1797            u=u0
1798            p=p0
1799            #
1800            #   initialize things:
1801            #
1802            converged=False
1803            #
1804            #  start loop:
1805            #
1806            #  initial search direction is g
1807            #
1808            while not converged :
1809                if self.iter>iter_max:
1810                    raise ArithmeticError("no convergence after %s steps."%self.iter)
1811                f_new=self.solve_f(u,p,tol)
1812                norm_f_new = util.Lsup(f_new)
1813                u_new=u-f_new
1814                g_new=self.solve_g(u_new,tol)
1815                self.iter+=1
1816                norm_g_new = util.sqrt(self.inner(g_new,g_new))
1817                if norm_f_new==0. and norm_g_new==0.: return u, p
1818                if self.iter>1 and not accepted_reduction==None:
1819                   #
1820                   #   did we manage to reduce the norm of G? I
1821                   #   if not we start a backtracking procedure
1822                   #
1823                   # print "new/old norm = ",norm_g_new, norm_g, norm_g_new/norm_g
1824                   if norm_g_new > accepted_reduction * norm_g:
1825                      sub_iter=0
1826                      s=self.relaxation
1827                      d=g
1828                      g_last=g
1829                      self.trace("    start substepping: f = %s, g = %s, relaxation = %s."%(norm_f_new, norm_g_new, s))
1830                      while sub_iter < self.subiter_max and  norm_g_new > accepted_reduction * norm_g:
1831                         dg= g_new-g_last
1832                         norm_dg=abs(util.sqrt(self.inner(dg,dg))/self.relaxation)
1833                         rad=self.inner(g_new,dg)/self.relaxation
1834                         # print "   ",sub_iter,": rad, norm_dg:",abs(rad), norm_dg*norm_g_new * angle_limit
1835                         # print "   ",sub_iter,": rad, norm_dg:",rad, norm_dg, norm_g_new, norm_g
1836                         if abs(rad) < norm_dg*norm_g_new * angle_limit:
1837                             if sub_iter>0: self.trace("    no further improvements expected from backtracking.")
1838                             break
1839                         r=self.relaxation
1840                         self.relaxation= - rad/norm_dg**2
1841                         s+=self.relaxation
1842                         #####
1843                         # a=g_new+self.relaxation*dg/r
1844                         # print "predicted new norm = ",util.sqrt(self.inner(a,a)),util.sqrt(self.inner(g_new,g_new)), self.relaxation
1845                         #####
1846                         g_last=g_new
1847                         p+=self.relaxation*d
1848                         f_new=self.solve_f(u,p,tol)
1849                         u_new=u-f_new
1850                         g_new=self.solve_g(u_new,tol)
1851                         self.iter+=1
1852                         norm_f_new = util.Lsup(f_new)
1853                         norm_g_new = util.sqrt(self.inner(g_new,g_new))
1854                         # print "   ",sub_iter," new g norm",norm_g_new
1855                         self.trace("    %s th sub-step: f = %s, g = %s, relaxation = %s."%(sub_iter, norm_f_new, norm_g_new, s))
1856                         #
1857                         #   can we expect reduction of g?
1858                         #
1859                         # u_last=u_new
1860                         sub_iter+=1
1861                      self.relaxation=s
1862                #
1863                #  check for convergence:
1864                #
1865                norm_u_new = util.Lsup(u_new)
1866                p_new=p+self.relaxation*g_new
1867                norm_p_new = util.sqrt(self.inner(p_new,p_new))
1868                self.trace("%s th step: f = %s, f/u = %s, g = %s, g/p = %s, relaxation = %s."%(self.iter, norm_f_new ,norm_f_new/norm_u_new, norm_g_new, norm_g_new/norm_p_new, self.relaxation))
1869    
1870                if self.iter>1:
1871                   dg2=g_new-g
1872                   df2=f_new-f
1873                   norm_dg2=util.sqrt(self.inner(dg2,dg2))
1874                   norm_df2=util.Lsup(df2)
1875                   # print norm_g_new, norm_g, norm_dg, norm_p, tolerance
1876                   tol_eq_g=tolerance*norm_dg2/(norm_g*abs(self.relaxation))*norm_p_new
1877                   tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new
1878                   if norm_g_new <= tol_eq_g and norm_f_new <= tol_eq_f:
1879                       converged=True
1880                f, norm_f, u, norm_u, g, norm_g, p, norm_p = f_new, norm_f_new, u_new, norm_u_new, g_new, norm_g_new, p_new, norm_p_new
1881            self.trace("convergence after %s steps."%self.iter)
1882            return u,p

Legend:
Removed from v.782  
changed lines
  Added in v.2100

  ViewVC Help
Powered by ViewVC 1.1.26