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

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

  ViewVC Help
Powered by ViewVC 1.1.26