/[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 1312 by ksteube, Mon Sep 24 06:18:44 2007 UTC revision 1878 by gross, Tue Oct 14 03:39:13 2008 UTC
# Line 1  Line 1 
1    
2    ########################################################
3  #  #
4  # $Id$  # Copyright (c) 2003-2008 by University of Queensland
5  #  # Earth Systems Science Computational Center (ESSCC)
6  #######################################################  # http://www.uq.edu.au/esscc
 #  
 #           Copyright 2003-2007 by ACceSS MNRF  
 #       Copyright 2007 by University of Queensland  
 #  
 #                http://esscc.uq.edu.au  
 #        Primary Business: Queensland, Australia  
 #  Licensed under the Open Software License version 3.0  
 #     http://www.opensource.org/licenses/osl-3.0.php  
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 32  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
# Line 48  import numarray Line 45  import numarray
45  import util  import util
46  import math  import math
47    
48    ##### Added by Artak
49    # from Numeric import zeros,Int,Float64
50    ###################################
51    
52    
53  class TimeIntegrationManager:  class TimeIntegrationManager:
54    """    """
55    a simple mechanism to manage time dependend values.    a simple mechanism to manage time dependend values.
# Line 435  class NegativeNorm(SolverSchemeException Line 437  class NegativeNorm(SolverSchemeException
437     """     """
438     pass     pass
439    
440  def PCG(b,x,Aprod,Msolve,bilinearform, norm, verbose=True, iter_max=100, tolerance=math.sqrt(util.EPSILON)):  class IterationHistory(object):
441     """     """
442     Solver for     The IterationHistory class is used to define a stopping criterium. It keeps track of the
443       residual norms. The stoppingcriterium indicates termination if the residual norm has been reduced by
444       a given tolerance.
445       """
446       def __init__(self,tolerance=math.sqrt(util.EPSILON),verbose=False):
447          """
448          Initialization
449    
450     M{Ax=b}        @param tolerance: tolerance
451          @type tolerance: positive C{float}
452          @param verbose: switches on the printing out some information
453          @type verbose: C{bool}
454          """
455          if not tolerance>0.:
456              raise ValueError,"tolerance needs to be positive."
457          self.tolerance=tolerance
458          self.verbose=verbose
459          self.history=[]
460       def stoppingcriterium(self,norm_r,r,x):
461           """
462           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.
463    
464     with a symmetric and positive definite operator A (more details required!).        
465     It uses the conjugate gradient method with preconditioner M providing an approximation of A.         @param norm_r: current residual norm
466           @type norm_r: non-negative C{float}
467           @param r: current residual (not used)
468           @param x: current solution approximation (not used)
469           @return: C{True} is the stopping criterium is fullfilled. Otherwise C{False} is returned.
470           @rtype: C{bool}
471    
472           """
473           self.history.append(norm_r)
474           if self.verbose: print "iter: %s:  inner(rhat,r) = %e"#(len(self.history)-1, self.history[-1])
475           return self.history[-1]<=self.tolerance * self.history[0]
476    
477       def stoppingcriterium2(self,norm_r,norm_b,solver="GMRES",TOL=None):
478           """
479           returns True if the C{norm_r} is C{tolerance}*C{norm_b}
480    
481          
482           @param norm_r: current residual norm
483           @type norm_r: non-negative C{float}
484           @param norm_b: norm of right hand side
485           @type norm_b: non-negative C{float}
486           @return: C{True} is the stopping criterium is fullfilled. Otherwise C{False} is returned.
487           @rtype: C{bool}
488    
489     The iteration is terminated if         """
490           if TOL==None:
491              TOL=self.tolerance
492           self.history.append(norm_r)
493           if self.verbose: print "iter: %s:  norm(r) = %e"#(len(self.history)-1, self.history[-1])
494           return self.history[-1]<=TOL * norm_b
495    
496     M{norm(r) <= tolerance * norm(b)}  def PCG(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
497       """
498       Solver for
499    
500     where C{norm()} defines a norm and     M{Ax=b}
501    
502     M{r = b-Ax}     with a symmetric and positive definite operator A (more details required!).
503       It uses the conjugate gradient method with preconditioner M providing an approximation of A.
504    
505     the residual.     The iteration is terminated if the C{stoppingcriterium} function return C{True}.
506    
507     For details on the preconditioned conjugate gradient method see the book:     For details on the preconditioned conjugate gradient method see the book:
508    
# Line 461  def PCG(b,x,Aprod,Msolve,bilinearform, n Line 511  def PCG(b,x,Aprod,Msolve,bilinearform, n
511     C. Romine, and H. van der Vorst.     C. Romine, and H. van der Vorst.
512    
513     @param b: the right hand side of the liner system. C{b} is altered.     @param b: the right hand side of the liner system. C{b} is altered.
514     @type b: any object type R supporting inplace add (x+=y) and scaling (x=scalar*y)     @type b: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
    @param x: an initial guess for the solution  
    @type x: any object type S supporting inplace add (x+=y), scaling (x=scalar*y)  
515     @param Aprod: returns the value Ax     @param Aprod: returns the value Ax
516     @type Aprod: function C{Aprod(x)} where C{x} is of object type S. The returned object needs to be of type R.     @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}.
517     @param Msolve: solves Mx=r     @param Msolve: solves Mx=r
518     @type Msolve: function C{Msolve(r)} where C{r} is of object type R. The returned object needs to be of tupe S.     @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
519    type like argument C{x}.
520     @param bilinearform: inner product C{<x,r>}     @param bilinearform: inner product C{<x,r>}
521     @type bilinearform: function C{bilinearform(x,r)} where C{x} is of object type S and C{r} is of object type R. The returned value is a C{float}.     @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}.
522     @param norm: norm calculation for the residual C{r=b-Ax}.     @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.
523     @type norm: function C{norm(r)} where C{r} is of object type R. The returned value is a C{float}.     @type stoppingcriterium: function that returns C{True} or C{False}
524     @param verbose: switches on the printing out some information     @param x: an initial guess for the solution. If no C{x} is given 0*b is used.
525     @type verbose: C{bool}     @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
526     @param iter_max: maximum number of iteration steps.     @param iter_max: maximum number of iteration steps.
527     @type iter_max: C{int}     @type iter_max: C{int}
528     @param tolerance: tolerance     @return: the solution approximation and the corresponding residual
529     @type tolerance: positive C{float}     @rtype: C{tuple}
530     @return: the solution apprximation and the corresponding residual     @warning: C{b} and C{x} are altered.
    @rtype: C{tuple} of an S type and and an R type object.A  
    @warning: C{b} ans C{x} are altered.  
531     """     """
    if verbose:  
         print "Enter PCG for solving Ax=b\n\t iter_max =%s\t tolerance   =%e"%(iter_max,tolerance)  
532     iter=0     iter=0
533     normb = norm(b)     if x==None:
534     if normb<0: raise NegativeNorm        x=0*b
535       else:
536     b += (-1)*Aprod(x)        b += (-1)*Aprod(x)
537     r=b     r=b
538     rhat=Msolve(r)     rhat=Msolve(r)
539     d = rhat;     d = rhat
540     rhat_dot_r = bilinearform(rhat, r)     rhat_dot_r = bilinearform(rhat, r)
541       if rhat_dot_r<0: raise NegativeNorm,"negative norm."
542    
543     while True:     while not stoppingcriterium(math.sqrt(rhat_dot_r),r,x):
544         normr=norm(r)         iter+=1
        if normr<0: raise NegativeNorm  
        if verbose: print "iter: %s: norm(r) = %e, tolerance*norm(b) = %e"%(iter, normr,tolerance * normb)  
        if normr <= tolerance * normb: return x,r  
   
        iter+=1 # k = iter = 1 first time through  
545         if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max         if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
546    
547         q=Aprod(d)         q=Aprod(d)
# Line 515  def PCG(b,x,Aprod,Msolve,bilinearform, n Line 556  def PCG(b,x,Aprod,Msolve,bilinearform, n
556         d=rhat         d=rhat
557    
558         rhat_dot_r = rhat_dot_r_new         rhat_dot_r = rhat_dot_r_new
559           if rhat_dot_r<0: raise NegativeNorm,"negative norm."
560    
561       return x,r
562    
563    class Defect(object):
564        """
565        defines a non-linear defect F(x) of a variable x
566        """
567        def __init__(self):
568            """
569            initialize defect
570            """
571            self.setDerivativeIncrementLength()
572    
573        def bilinearform(self, x0, x1):
574            """
575            returns the inner product of x0 and x1
576            @param x0: a value for x
577            @param x1: a value for x
578            @return: the inner product of x0 and x1
579            @rtype: C{float}
580            """
581            return 0
582          
583        def norm(self,x):
584            """
585            the norm of argument C{x}
586    
587            @param x: a value for x
588            @return: norm of argument x
589            @rtype: C{float}
590            @note: by default C{sqrt(self.bilinearform(x,x)} is retrurned.
591            """
592            s=self.bilinearform(x,x)
593            if s<0: raise NegativeNorm,"negative norm."
594            return math.sqrt(s)
595    
596    
597        def eval(self,x):
598            """
599            returns the value F of a given x
600    
601            @param x: value for which the defect C{F} is evalulated.
602            @return: value of the defect at C{x}
603            """
604            return 0
605    
606        def __call__(self,x):
607            return self.eval(x)
608    
609        def setDerivativeIncrementLength(self,inc=math.sqrt(util.EPSILON)):
610            """
611            sets the relative length of the increment used to approximate the derivative of the defect
612            the increment is inc*norm(x)/norm(v)*v in the direction of v with x as a staring point.
613    
614            @param inc: relative increment length
615            @type inc: positive C{float}
616            """
617            if inc<=0: raise ValueError,"positive increment required."
618            self.__inc=inc
619    
620        def getDerivativeIncrementLength(self):
621            """
622            returns the relative increment length used to approximate the derivative of the defect
623            @return: value of the defect at C{x}
624            @rtype: positive C{float}
625            """
626            return self.__inc
627    
628        def derivative(self, F0, x0, v, v_is_normalised=True):
629            """
630            returns the directional derivative at x0 in the direction of v
631    
632            @param F0: value of this defect at x0
633            @param x0: value at which derivative is calculated.
634            @param v: direction
635            @param v_is_normalised: is true to indicate that C{v} is nomalized (self.norm(v)=0)
636            @return: derivative of this defect at x0 in the direction of C{v}
637            @note: by default numerical evaluation (self.eval(x0+eps*v)-F0)/eps is used but this method
638            maybe oepsnew verwritten to use exact evalution.
639            """
640            normx=self.norm(x0)
641            if normx>0:
642                 epsnew = self.getDerivativeIncrementLength() * normx
643            else:
644                 epsnew = self.getDerivativeIncrementLength()
645            if not v_is_normalised:
646                normv=self.norm(v)
647                if normv<=0:
648                   return F0*0
649                else:
650                   epsnew /= normv
651            F1=self.eval(x0 + epsnew * v)
652            return (F1-F0)/epsnew
653    
654    ######################################    
655    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):
656       """
657       solves a non-linear problem M{F(x)=0} for unknown M{x} using the stopping criterion:
658    
659       M{norm(F(x) <= atol + rtol * norm(F(x0)}
660      
661       where M{x0} is the initial guess.
662    
663       @param defect: object defining the the function M{F}, C{defect.norm} defines the M{norm} used in the stopping criterion.
664       @type defect: L{Defect}
665       @param x: initial guess for the solution, C{x} is altered.
666       @type x: any object type allowing basic operations such as  L{numarray.NumArray}, L{Data}
667       @param iter_max: maximum number of iteration steps
668       @type iter_max: positive C{int}
669       @param sub_iter_max:
670       @type sub_iter_max:
671       @param atol: absolute tolerance for the solution
672       @type atol: positive C{float}
673       @param rtol: relative tolerance for the solution
674       @type rtol: positive C{float}
675       @param gamma: tolerance safety factor for inner iteration
676       @type gamma: positive C{float}, less than 1
677       @param sub_tol_max: upper bound for inner tolerance.
678       @type sub_tol_max: positive C{float}, less than 1
679       @return: an approximation of the solution with the desired accuracy
680       @rtype: same type as the initial guess.
681       """
682       lmaxit=iter_max
683       if atol<0: raise ValueError,"atol needs to be non-negative."
684       if rtol<0: raise ValueError,"rtol needs to be non-negative."
685       if rtol+atol<=0: raise ValueError,"rtol or atol needs to be non-negative."
686       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
687       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
688    
689       F=defect(x)
690       fnrm=defect.norm(F)
691       stop_tol=atol + rtol*fnrm
692       sub_tol=sub_tol_max
693       if verbose: print "NewtonGMRES: initial residual = %e."%fnrm
694       if verbose: print "             tolerance = %e."%sub_tol
695       iter=1
696       #
697       # main iteration loop
698       #
699       while not fnrm<=stop_tol:
700                if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
701                #
702            #   adjust sub_tol_
703            #
704                if iter > 1:
705               rat=fnrm/fnrmo
706                   sub_tol_old=sub_tol
707               sub_tol=gamma*rat**2
708               if gamma*sub_tol_old**2 > .1: sub_tol=max(sub_tol,gamma*sub_tol_old**2)
709               sub_tol=max(min(sub_tol,sub_tol_max), .5*stop_tol/fnrm)
710            #
711            # calculate newton increment xc
712                #     if iter_max in __FDGMRES is reached MaxIterReached is thrown
713                #     if iter_restart -1 is returned as sub_iter
714                #     if  atol is reached sub_iter returns the numer of steps performed to get there
715                #
716                #  
717                if verbose: print "             subiteration (GMRES) is called with relative tolerance %e."%sub_tol
718                try:
719                   xc, sub_iter=__FDGMRES(F, defect, x, sub_tol*fnrm, iter_max=iter_max-iter, iter_restart=sub_iter_max)
720                except MaxIterReached:
721                   raise MaxIterReached,"maximum number of %s steps reached."%iter_max
722                if sub_iter<0:
723                   iter+=sub_iter_max
724                else:
725                   iter+=sub_iter
726                # ====
727            x+=xc
728                F=defect(x)
729            iter+=1
730                fnrmo, fnrm=fnrm, defect.norm(F)
731                if verbose: print "             step %s: residual %e."%(iter,fnrm)
732       if verbose: print "NewtonGMRES: completed after %s steps."%iter
733       return x
734    
735    def __givapp(c,s,vin):
736        """
737        apply a sequence of Givens rotations (c,s) to the recuirsively to the vector vin
738        @warning: C{vin} is altered.
739        """
740        vrot=vin
741        if isinstance(c,float):
742            vrot=[c*vrot[0]-s*vrot[1],s*vrot[0]+c*vrot[1]]
743        else:
744            for i in range(len(c)):
745                w1=c[i]*vrot[i]-s[i]*vrot[i+1]
746            w2=s[i]*vrot[i]+c[i]*vrot[i+1]
747                vrot[i:i+2]=w1,w2
748        return vrot
749    
750    def __FDGMRES(F0, defect, x0, atol, iter_max=100, iter_restart=20):
751       h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)
752       c=numarray.zeros(iter_restart,numarray.Float64)
753       s=numarray.zeros(iter_restart,numarray.Float64)
754       g=numarray.zeros(iter_restart,numarray.Float64)
755       v=[]
756    
757       rho=defect.norm(F0)
758       if rho<=0.: return x0*0
759      
760       v.append(-F0/rho)
761       g[0]=rho
762       iter=0
763       while rho > atol and iter<iter_restart-1:
764    
765        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
766    
767            p=defect.derivative(F0,x0,v[iter], v_is_normalised=True)
768        v.append(p)
769    
770        v_norm1=defect.norm(v[iter+1])
771    
772            # Modified Gram-Schmidt
773        for j in range(iter+1):
774             h[j][iter]=defect.bilinearform(v[j],v[iter+1])  
775             v[iter+1]-=h[j][iter]*v[j]
776          
777        h[iter+1][iter]=defect.norm(v[iter+1])
778        v_norm2=h[iter+1][iter]
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=defect.bilinearform(v[j],v[iter+1])
784                   h[j][iter]=h[j][iter]+hr
785                   v[iter+1] -= hr*v[j]
786    
787            v_norm2=defect.norm(v[iter+1])
788            h[iter+1][iter]=v_norm2
789            #   watch out for happy breakdown
790            if not v_norm2 == 0:
791                    v[iter+1]=v[iter+1]/h[iter+1][iter]
792    
793            #   Form and store the information for the new Givens rotation
794        if iter > 0 :
795            hhat=numarray.zeros(iter+1,numarray.Float64)
796            for i in range(iter+1) : hhat[i]=h[i][iter]
797            hhat=__givapp(c[0:iter],s[0:iter],hhat);
798                for i in range(iter+1) : h[i][iter]=hhat[i]
799    
800        mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])
801    
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       # At this point either iter > iter_max or rho < tol.
814       # It's time to compute x and leave.        
815       if iter > 0 :
816         y=numarray.zeros(iter,numarray.Float64)    
817         y[iter-1] = g[iter-1] / h[iter-1][iter-1]
818         if iter > 1 :  
819            i=iter-2  
820            while i>=0 :
821              y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]
822              i=i-1
823         xhat=v[iter-1]*y[iter-1]
824         for i in range(iter-1):
825        xhat += v[i]*y[i]
826       else :
827          xhat=v[0] * 0
828    
829       if iter<iter_restart-1:
830          stopped=iter
831       else:
832          stopped=-1
833    
834       return xhat,stopped
835    
836    ##############################################
837    def GMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):
838    ################################################
839       m=iter_restart
840       iter=0
841       xc=x
842       while True:
843          if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached"%iter_max
844          xc,stopped=__GMRESm(b*1, Aprod, Msolve, bilinearform, stoppingcriterium, x=xc*1, iter_max=iter_max-iter, iter_restart=m)
845          if stopped: break
846          iter+=iter_restart    
847       return xc
848    
849    def __GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):
850       iter=0
851       r=Msolve(b)
852       r_dot_r = bilinearform(r, r)
853       if r_dot_r<0: raise NegativeNorm,"negative norm."
854       norm_b=math.sqrt(r_dot_r)
855    
856       if x==None:
857          x=0*b
858       else:
859          r=Msolve(b-Aprod(x))
860          r_dot_r = bilinearform(r, r)
861          if r_dot_r<0: raise NegativeNorm,"negative norm."
862      
863       h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)
864       c=numarray.zeros(iter_restart,numarray.Float64)
865       s=numarray.zeros(iter_restart,numarray.Float64)
866       g=numarray.zeros(iter_restart,numarray.Float64)
867       v=[]
868    
869       rho=math.sqrt(r_dot_r)
870      
871       v.append(r/rho)
872       g[0]=rho
873    
874       while not (stoppingcriterium(rho,norm_b) or iter==iter_restart-1):
875    
876        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
877    
878        p=Msolve(Aprod(v[iter]))
879    
880        v.append(p)
881    
882        v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
883    
884    # Modified Gram-Schmidt
885        for j in range(iter+1):
886          h[j][iter]=bilinearform(v[j],v[iter+1])  
887          v[iter+1]-=h[j][iter]*v[j]
888          
889        h[iter+1][iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
890        v_norm2=h[iter+1][iter]
891    
892    # Reorthogonalize if needed
893        if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
894         for j in range(iter+1):  
895            hr=bilinearform(v[j],v[iter+1])
896                h[j][iter]=h[j][iter]+hr
897                v[iter+1] -= hr*v[j]
898    
899         v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
900         h[iter+1][iter]=v_norm2
901    
902    #   watch out for happy breakdown
903            if not v_norm2 == 0:
904             v[iter+1]=v[iter+1]/h[iter+1][iter]
905    
906    #   Form and store the information for the new Givens rotation
907        if iter > 0 :
908            hhat=numarray.zeros(iter+1,numarray.Float64)
909            for i in range(iter+1) : hhat[i]=h[i][iter]
910            hhat=__givapp(c[0:iter],s[0:iter],hhat);
911                for i in range(iter+1) : h[i][iter]=hhat[i]
912    
913        mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])
914    
915        if mu!=0 :
916            c[iter]=h[iter][iter]/mu
917            s[iter]=-h[iter+1][iter]/mu
918            h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter]
919            h[iter+1][iter]=0.0
920            g[iter:iter+2]=__givapp(c[iter],s[iter],g[iter:iter+2])
921    
922    # Update the residual norm
923                  
924            rho=abs(g[iter+1])
925        iter+=1
926    
927    # At this point either iter > iter_max or rho < tol.
928    # It's time to compute x and leave.        
929    
930       if iter > 0 :
931         y=numarray.zeros(iter,numarray.Float64)    
932         y[iter-1] = g[iter-1] / h[iter-1][iter-1]
933         if iter > 1 :  
934            i=iter-2  
935            while i>=0 :
936              y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]
937              i=i-1
938         xhat=v[iter-1]*y[iter-1]
939         for i in range(iter-1):
940        xhat += v[i]*y[i]
941       else : xhat=v[0]
942    
943       x += xhat
944       if iter<iter_restart-1:
945          stopped=True
946       else:
947          stopped=False
948    
949       return x,stopped
950    
951    #################################################
952    def MINRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
953    #################################################
954        #
955        #  minres solves the system of linear equations Ax = b
956        #  where A is a symmetric matrix (possibly indefinite or singular)
957        #  and b is a given vector.
958        #  
959        #  "A" may be a dense or sparse matrix (preferably sparse!)
960        #  or the name of a function such that
961        #               y = A(x)
962        #  returns the product y = Ax for any given vector x.
963        #
964        #  "M" defines a positive-definite preconditioner M = C C'.
965        #  "M" may be a dense or sparse matrix (preferably sparse!)
966        #  or the name of a function such that
967        #  solves the system My = x for any given vector x.
968        #
969        #
970        
971        #------------------------------------------------------------------
972        # Set up y and v for the first Lanczos vector v1.
973        # y  =  beta1 P' v1,  where  P = C**(-1).
974        # v is really P' v1.
975        #------------------------------------------------------------------
976        if x==None:
977          x=0*b
978        else:
979          b += (-1)*Aprod(x)
980    
981        r1    = b
982        y = Msolve(b)
983        beta1 = bilinearform(y,b)
984    
985        if beta1< 0: raise NegativeNorm,"negative norm."
986    
987        #  If b = 0 exactly, stop with x = 0.
988        if beta1==0: return x*0.
989    
990        if beta1> 0:
991          beta1  = math.sqrt(beta1)      
992    
993        #------------------------------------------------------------------
994        # Initialize quantities.
995        # ------------------------------------------------------------------
996        iter   = 0
997        Anorm = 0
998        ynorm = 0
999        oldb   = 0
1000        beta   = beta1
1001        dbar   = 0
1002        epsln  = 0
1003        phibar = beta1
1004        rhs1   = beta1
1005        rhs2   = 0
1006        rnorm  = phibar
1007        tnorm2 = 0
1008        ynorm2 = 0
1009        cs     = -1
1010        sn     = 0
1011        w      = b*0.
1012        w2     = b*0.
1013        r2     = r1
1014        eps    = 0.0001
1015    
1016        #---------------------------------------------------------------------
1017        # Main iteration loop.
1018        # --------------------------------------------------------------------
1019        while not stoppingcriterium(rnorm,Anorm*ynorm,'MINRES'):    #  checks ||r|| < (||A|| ||x||) * TOL
1020    
1021        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
1022            iter    = iter  +  1
1023    
1024            #-----------------------------------------------------------------
1025            # Obtain quantities for the next Lanczos vector vk+1, k = 1, 2,...
1026            # The general iteration is similar to the case k = 1 with v0 = 0:
1027            #
1028            #   p1      = Operator * v1  -  beta1 * v0,
1029            #   alpha1  = v1'p1,
1030            #   q2      = p2  -  alpha1 * v1,
1031            #   beta2^2 = q2'q2,
1032            #   v2      = (1/beta2) q2.
1033            #
1034            # Again, y = betak P vk,  where  P = C**(-1).
1035            #-----------------------------------------------------------------
1036            s = 1/beta                 # Normalize previous vector (in y).
1037            v = s*y                    # v = vk if P = I
1038        
1039            y      = Aprod(v)
1040        
1041            if iter >= 2:
1042              y = y - (beta/oldb)*r1
1043    
1044            alfa   = bilinearform(v,y)              # alphak
1045            y      += (- alfa/beta)*r2
1046            r1     = r2
1047            r2     = y
1048            y = Msolve(r2)
1049            oldb   = beta                           # oldb = betak
1050            beta   = bilinearform(y,r2)             # beta = betak+1^2
1051            if beta < 0: raise NegativeNorm,"negative norm."
1052    
1053            beta   = math.sqrt( beta )
1054            tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta
1055            
1056            if iter==1:                 # Initialize a few things.
1057              gmax   = abs( alfa )      # alpha1
1058              gmin   = gmax             # alpha1
1059    
1060            # Apply previous rotation Qk-1 to get
1061            #   [deltak epslnk+1] = [cs  sn][dbark    0   ]
1062            #   [gbar k dbar k+1]   [sn -cs][alfak betak+1].
1063        
1064            oldeps = epsln
1065            delta  = cs * dbar  +  sn * alfa  # delta1 = 0         deltak
1066            gbar   = sn * dbar  -  cs * alfa  # gbar 1 = alfa1     gbar k
1067            epsln  =               sn * beta  # epsln2 = 0         epslnk+1
1068            dbar   =            -  cs * beta  # dbar 2 = beta2     dbar k+1
1069    
1070            # Compute the next plane rotation Qk
1071    
1072            gamma  = math.sqrt(gbar*gbar+beta*beta)  # gammak
1073            gamma  = max(gamma,eps)
1074            cs     = gbar / gamma             # ck
1075            sn     = beta / gamma             # sk
1076            phi    = cs * phibar              # phik
1077            phibar = sn * phibar              # phibark+1
1078    
1079            # Update  x.
1080    
1081            denom = 1/gamma
1082            w1    = w2
1083            w2    = w
1084            w     = (v - oldeps*w1 - delta*w2) * denom
1085            x     +=  phi*w
1086    
1087            # Go round again.
1088    
1089            gmax   = max(gmax,gamma)
1090            gmin   = min(gmin,gamma)
1091            z      = rhs1 / gamma
1092            ynorm2 = z*z  +  ynorm2
1093            rhs1   = rhs2 -  delta*z
1094            rhs2   =      -  epsln*z
1095    
1096            # Estimate various norms and test for convergence.
1097    
1098            Anorm  = math.sqrt( tnorm2 )
1099            ynorm  = math.sqrt( ynorm2 )
1100    
1101            rnorm  = phibar
1102    
1103        return x
1104    
1105    def TFQMR(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
1106    
1107    # TFQMR solver for linear systems
1108    #
1109    #
1110    # initialization
1111    #
1112      errtol = math.sqrt(bilinearform(b,b))
1113      norm_b=errtol
1114      kmax  = iter_max
1115      error = []
1116    
1117      if math.sqrt(bilinearform(x,x)) != 0.0:
1118        r = b - Aprod(x)
1119      else:
1120        r = b
1121    
1122      r=Msolve(r)
1123    
1124      u1=0
1125      u2=0
1126      y1=0
1127      y2=0
1128    
1129      w = r
1130      y1 = r
1131      iter = 0
1132      d = 0
1133      
1134      v = Msolve(Aprod(y1))
1135      u1 = v
1136      
1137      theta = 0.0;
1138      eta = 0.0;
1139      tau = math.sqrt(bilinearform(r,r))
1140      error = [ error, tau ]
1141      rho = tau * tau
1142      m=1
1143    #
1144    #  TFQMR iteration
1145    #
1146    #  while ( iter < kmax-1 ):
1147      
1148      while not stoppingcriterium(tau*math.sqrt ( m + 1 ),norm_b,'TFQMR'):
1149        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
1150    
1151        sigma = bilinearform(r,v)
1152    
1153        if ( sigma == 0.0 ):
1154          raise 'TFQMR breakdown, sigma=0'
1155        
1156    
1157        alpha = rho / sigma
1158    
1159        for j in range(2):
1160    #
1161    #   Compute y2 and u2 only if you have to
1162    #
1163          if ( j == 1 ):
1164            y2 = y1 - alpha * v
1165            u2 = Msolve(Aprod(y2))
1166    
1167          m = 2 * (iter+1) - 2 + (j+1)
1168          if j==0:
1169             w = w - alpha * u1
1170             d = y1 + ( theta * theta * eta / alpha ) * d
1171          if j==1:
1172             w = w - alpha * u2
1173             d = y2 + ( theta * theta * eta / alpha ) * d
1174    
1175          theta = math.sqrt(bilinearform(w,w))/ tau
1176          c = 1.0 / math.sqrt ( 1.0 + theta * theta )
1177          tau = tau * theta * c
1178          eta = c * c * alpha
1179          x = x + eta * d
1180    #
1181    #  Try to terminate the iteration at each pass through the loop
1182    #
1183         # if ( tau * math.sqrt ( m + 1 ) <= errtol ):
1184         #   error = [ error, tau ]
1185         #   total_iters = iter
1186         #   break
1187          
1188    
1189        if ( rho == 0.0 ):
1190          raise 'TFQMR breakdown, rho=0'
1191        
1192    
1193        rhon = bilinearform(r,w)
1194        beta = rhon / rho;
1195        rho = rhon;
1196        y1 = w + beta * y2;
1197        u1 = Msolve(Aprod(y1))
1198        v = u1 + beta * ( u2 + beta * v )
1199        error = [ error, tau ]
1200        total_iters = iter
1201        
1202        iter = iter + 1
1203    
1204      return x
1205    
1206    
1207    #############################################
1208    
1209    class ArithmeticTuple(object):
1210       """
1211       tuple supporting inplace update x+=y and scaling x=a*y where x,y is an ArithmeticTuple and a is a float.
1212    
1213       example of usage:
1214    
1215       from esys.escript import Data
1216       from numarray import array
1217       a=Data(...)
1218       b=array([1.,4.])
1219       x=ArithmeticTuple(a,b)
1220       y=5.*x
1221    
1222       """
1223       def __init__(self,*args):
1224           """
1225           initialize object with elements args.
1226    
1227           @param args: tuple of object that support implace add (x+=y) and scaling (x=a*y)
1228           """
1229           self.__items=list(args)
1230    
1231       def __len__(self):
1232           """
1233           number of items
1234    
1235           @return: number of items
1236           @rtype: C{int}
1237           """
1238           return len(self.__items)
1239    
1240       def __getitem__(self,index):
1241           """
1242           get an item
1243    
1244           @param index: item to be returned
1245           @type index: C{int}
1246           @return: item with index C{index}
1247           """
1248           return self.__items.__getitem__(index)
1249    
1250       def __mul__(self,other):
1251           """
1252           scaling from the right
1253    
1254           @param other: scaling factor
1255           @type other: C{float}
1256           @return: itemwise self*other
1257           @rtype: L{ArithmeticTuple}
1258           """
1259           out=[]
1260           other=1.*other
1261           if isinstance(other,float):
1262        for i in range(len(self)):
1263               out.append(self[i]*other)
1264           else:
1265            for i in range(len(self)):
1266               out.append(self[i]*other[i])
1267           return ArithmeticTuple(*tuple(out))
1268    
1269       def __rmul__(self,other):
1270           """
1271           scaling from the left
1272    
1273           @param other: scaling factor
1274           @type other: C{float}
1275           @return: itemwise other*self
1276           @rtype: L{ArithmeticTuple}
1277           """
1278           out=[]
1279           other=1.*other
1280           if isinstance(other,float):
1281        for i in range(len(self)):
1282               out.append(other*self[i])
1283           else:
1284            for i in range(len(self)):
1285               out.append(other[i]*self[i])
1286           return ArithmeticTuple(*tuple(out))
1287    
1288    #########################
1289    # Added by Artak
1290    #########################
1291       def __div__(self,other):
1292           """
1293           dividing from the right
1294    
1295           @param other: scaling factor
1296           @type other: C{float}
1297           @return: itemwise self/other
1298           @rtype: L{ArithmeticTuple}
1299           """
1300           out=[]
1301           other=1.*other
1302           if isinstance(other,float):
1303        for i in range(len(self)):
1304               out.append(self[i]*(1./other))
1305           else:
1306            for i in range(len(self)):
1307               out.append(self[i]*(1./other[i]))
1308           return ArithmeticTuple(*tuple(out))
1309    
1310       def __rdiv__(self,other):
1311           """
1312           dividing from the left
1313    
1314           @param other: scaling factor
1315           @type other: C{float}
1316           @return: itemwise other/self
1317           @rtype: L{ArithmeticTuple}
1318           """
1319           out=[]
1320           other=1.*other
1321           if isinstance(other,float):
1322            for i in range(len(self)):
1323               out.append(other*(1./self[i]))
1324           else:
1325            for i in range(len(self)):
1326               out.append(other[i]*(1./self[i]))
1327           return ArithmeticTuple(*tuple(out))
1328      
1329    ##########################################33
1330    
1331       def __iadd__(self,other):
1332           """
1333           in-place add of other to self
1334    
1335           @param other: increment
1336           @type other: C{ArithmeticTuple}
1337           """
1338           if len(self) != len(other):
1339               raise ValueError,"tuple length must match."
1340           for i in range(len(self)):
1341               self.__items[i]+=other[i]
1342           return self
1343    
1344       def __add__(self,other):
1345           """
1346           add other to self
1347    
1348           @param other: increment
1349           @type other: C{ArithmeticTuple}
1350           """
1351    #      if not isinstance(other,float):
1352           if len(self) != len(other):
1353              raise ValueError,"tuple length must match."
1354           for i in range(len(self)):
1355              self.__items[i]+=other[i]
1356    #       else:
1357    #        for i in range(len(self)):
1358    #           self.__items[i]+=other
1359    
1360           return self
1361    
1362       def __sub__(self,other):
1363           """
1364           subtract other from self
1365    
1366           @param other: increment
1367           @type other: C{ArithmeticTuple}
1368           """
1369           if len(self) != len(other):
1370               raise ValueError,"tuple length must match."
1371           for i in range(len(self)):
1372               self.__items[i]-=other[i]
1373           return self
1374      
1375       def __isub__(self,other):
1376           """
1377           subtract other from self
1378    
1379           @param other: increment
1380           @type other: C{ArithmeticTuple}
1381           """
1382           if len(self) != len(other):
1383               raise ValueError,"tuple length must match."
1384           for i in range(len(self)):
1385               self.__items[i]-=other[i]
1386           return self
1387    
1388       def __neg__(self):
1389           """
1390           negate
1391    
1392           """
1393           for i in range(len(self)):
1394               self.__items[i]=-self.__items[i]
1395           return self
1396    
1397    
1398    class HomogeneousSaddlePointProblem(object):
1399          """
1400          This provides a framwork for solving linear homogeneous saddle point problem of the form
1401    
1402                 Av+B^*p=f
1403                 Bv    =0
1404    
1405          for the unknowns v and p and given operators A and B and given right hand side f.
1406          B^* is the adjoint operator of B is the given inner product.
1407    
1408          """
1409          def __init__(self,**kwargs):
1410            self.setTolerance()
1411            self.setToleranceReductionFactor()
1412    
1413          def initialize(self):
1414            """
1415            initialize the problem (overwrite)
1416            """
1417            pass
1418          def B(self,v):
1419             """
1420             returns Bv (overwrite)
1421             @rtype: equal to the type of p
1422    
1423             @note: boundary conditions on p should be zero!
1424             """
1425             pass
1426    
1427          def inner(self,p0,p1):
1428             """
1429             returns inner product of two element p0 and p1  (overwrite)
1430            
1431             @type p0: equal to the type of p
1432             @type p1: equal to the type of p
1433             @rtype: C{float}
1434    
1435             @rtype: equal to the type of p
1436             """
1437             pass
1438    
1439          def solve_A(self,u,p):
1440             """
1441             solves Av=f-Au-B^*p with accuracy self.getReducedTolerance() (overwrite)
1442    
1443             @rtype: equal to the type of v
1444             @note: boundary conditions on v should be zero!
1445             """
1446             pass
1447    
1448          def solve_prec(self,p):
1449             """
1450             provides a preconditioner for BA^{-1}B^* with accuracy self.getReducedTolerance() (overwrite)
1451    
1452             @rtype: equal to the type of p
1453             """
1454             pass
1455    
1456          def stoppingcriterium(self,Bv,v,p):
1457             """
1458             returns a True if iteration is terminated. (overwrite)
1459    
1460             @rtype: C{bool}
1461             """
1462             pass
1463                
1464          def __inner(self,p,r):
1465             return self.inner(p,r[1])
1466    
1467          def __inner_p(self,p1,p2):
1468             return self.inner(p1,p2)
1469          
1470          def __inner_a(self,a1,a2):
1471             return self.inner_a(a1,a2)
1472    
1473          def __inner_a1(self,a1,a2):
1474             return self.inner(a1[1],a2[1])
1475    
1476          def __stoppingcriterium(self,norm_r,r,p):
1477              return self.stoppingcriterium(r[1],r[0],p)
1478    
1479          def __stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None):
1480              return self.stoppingcriterium2(norm_r,norm_b,solver,TOL)
1481    
1482          def setTolerance(self,tolerance=1.e-8):
1483                  self.__tol=tolerance
1484          def getTolerance(self):
1485                  return self.__tol
1486          def setToleranceReductionFactor(self,reduction=0.01):
1487                  self.__reduction=reduction
1488          def getSubProblemTolerance(self):
1489                  return self.__reduction*self.getTolerance()
1490    
1491          def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='PCG',iter_restart=20):
1492                  """
1493                  solves the saddle point problem using initial guesses v and p.
1494    
1495                  @param max_iter: maximum number of iteration steps.
1496                  """
1497                  self.verbose=verbose
1498                  self.show_details=show_details and self.verbose
1499    
1500                  # assume p is known: then v=A^-1(f-B^*p)
1501                  # which leads to BA^-1B^*p = BA^-1f  
1502    
1503              # Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask)      
1504              self.__z=v+self.solve_A(v,p*0)
1505                  Bz=self.B(self.__z)
1506                  #
1507              #   solve BA^-1B^*p = Bz
1508                  #
1509                  #
1510                  #
1511                  self.iter=0
1512              if solver=='GMRES':      
1513                    if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter
1514                    p=GMRES(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.,iter_restart=iter_restart)
1515                    # solve Au=f-B^*p
1516                    #       A(u-v)=f-B^*p-Av
1517                    #       u=v+(u-v)
1518            u=v+self.solve_A(v,p)
1519    
1520              if solver=='TFQMR':      
1521                    if self.verbose: print "enter TFQMR method (iter_max=%s)"%max_iter
1522                    p=TFQMR(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)
1523                    # solve Au=f-B^*p
1524                    #       A(u-v)=f-B^*p-Av
1525                    #       u=v+(u-v)
1526            u=v+self.solve_A(v,p)
1527    
1528              if solver=='MINRES':      
1529                    if self.verbose: print "enter MINRES method (iter_max=%s)"%max_iter
1530                    p=MINRES(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)
1531                    # solve Au=f-B^*p
1532                    #       A(u-v)=f-B^*p-Av
1533                    #       u=v+(u-v)
1534            u=v+self.solve_A(v,p)
1535                  
1536              if solver=='GMRESC':      
1537                    if self.verbose: print "enter GMRES coupled method (iter_max=%s)"%max_iter
1538                    p0=self.solve_prec1(Bz)
1539                #alfa=(1/self.vol)*util.integrate(util.interpolate(p,escript.Function(self.domain)))
1540                    #p-=alfa
1541                    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)
1542                    #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())
1543    
1544                    # solve Au=f-B^*p
1545                    #       A(u-v)=f-B^*p-Av
1546                    #       u=v+(u-v)
1547                p=x[1]
1548            u=v+self.solve_A(v,p)      
1549            #p=x[1]
1550            #u=x[0]
1551    
1552                  if solver=='PCG':
1553                    #   note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv
1554                    #
1555                    #   with                    Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask)
1556                    #                           A(v-z)= f -Az - B^*p (v-z=0 on fixed_u_mask)
1557                    if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter
1558                    p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p)
1559                u=r[0]  
1560                    # print "DIFF=",util.integrate(p)
1561    
1562                  # print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)
1563    
1564              return u,p
1565    
1566          def __Msolve(self,r):
1567              return self.solve_prec(r[1])
1568    
1569          def __Msolve2(self,r):
1570              return self.solve_prec(r*1)
1571    
1572          def __Mempty(self,r):
1573              return r
1574    
1575    
1576          def __Aprod(self,p):
1577              # return BA^-1B*p
1578              #solve Av =B^*p as Av =f-Az-B^*(-p)
1579              v=self.solve_A(self.__z,-p)
1580              return ArithmeticTuple(v, self.B(v))
1581    
1582          def __Anext(self,x):
1583              # return next v,p
1584              #solve Av =-B^*p as Av =f-Az-B^*p
1585    
1586          pc=x[1]
1587              v=self.solve_A(self.__z,-pc)
1588          p=self.solve_prec1(self.B(v))
1589    
1590              return ArithmeticTuple(v,p)
1591    
1592    
1593          def __Aprod2(self,p):
1594              # return BA^-1B*p
1595              #solve Av =B^*p as Av =f-Az-B^*(-p)
1596          v=self.solve_A(self.__z,-p)
1597              return self.B(v)
1598    
1599          def __Aprod_Newton(self,p):
1600              # return BA^-1B*p - Bz
1601              #solve Av =-B^*p as Av =f-Az-B^*p
1602          v=self.solve_A(self.__z,-p)
1603              return self.B(v-self.__z)
1604    
1605          def __Aprod_Newton2(self,x):
1606              # return BA^-1B*p - Bz
1607              #solve Av =-B^*p as Av =f-Az-B^*p
1608              pc=x[1]
1609          v=self.solve_A(self.__z,-pc)
1610              p=self.solve_prec1(self.B(v-self.__z))
1611              return ArithmeticTuple(v,p)
1612    
1613    
1614    def MaskFromBoundaryTag(domain,*tags):
1615       """
1616       create a mask on the Solution(domain) function space which one for samples
1617       that touch the boundary tagged by tags.
1618    
1619       usage: m=MaskFromBoundaryTag(domain,"left", "right")
1620    
1621       @param domain: a given domain
1622       @type domain: L{escript.Domain}
1623       @param tags: boundray tags
1624       @type tags: C{str}
1625       @return: a mask which marks samples that are touching the boundary tagged by any of the given tags.
1626       @rtype: L{escript.Data} of rank 0
1627       """
1628       pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1629       d=escript.Scalar(0.,escript.FunctionOnBoundary(function_space.getDomain()))
1630       for t in tags: d.setTaggedValue(t,1.)
1631       pde.setValue(y=d)
1632       return util.whereNonZero(pde.getRightHandSide())
1633    #============================================================================================================================================
1634  class SaddlePointProblem(object):  class SaddlePointProblem(object):
1635     """     """
1636     This implements a solver for a saddlepoint problem     This implements a solver for a saddlepoint problem
# Line 542  class SaddlePointProblem(object): Line 1656  class SaddlePointProblem(object):
1656         @type verbose: C{bool}         @type verbose: C{bool}
1657         @note: this method may be overwritten by a particular saddle point problem         @note: this method may be overwritten by a particular saddle point problem
1658         """         """
1659           print "SaddlePointProblem should not be used anymore!"
1660         if not isinstance(verbose,bool):         if not isinstance(verbose,bool):
1661              raise TypeError("verbose needs to be of type bool.")              raise TypeError("verbose needs to be of type bool.")
1662         self.__verbose=verbose         self.__verbose=verbose
1663         self.relaxation=1.         self.relaxation=1.
1664           DeprecationWarning("SaddlePointProblem should not be used anymore.")
1665    
1666     def trace(self,text):     def trace(self,text):
1667         """         """
# Line 719  class SaddlePointProblem(object): Line 1835  class SaddlePointProblem(object):
1835              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              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
1836          self.trace("convergence after %s steps."%self.iter)          self.trace("convergence after %s steps."%self.iter)
1837          return u,p          return u,p
 #   def solve(self,u0,p0,tolerance=1.e-6,iter_max=10,self.relaxation=1.):  
 #      tol=1.e-2  
 #      iter=0  
 #      converged=False  
 #      u=u0*1.  
 #      p=p0*1.  
 #      while not converged and iter<iter_max:  
 #          du=self.solve_f(u,p,tol)  
 #          u-=du  
 #          norm_du=util.Lsup(du)  
 #          norm_u=util.Lsup(u)  
 #          
 #          dp=self.relaxation*self.solve_g(u,tol)  
 #          p+=dp  
 #          norm_dp=util.sqrt(self.inner(dp,dp))  
 #          norm_p=util.sqrt(self.inner(p,p))  
 #          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)  
 #          iter+=1  
 #  
 #          converged = (norm_du <= tolerance*norm_u) and  (norm_dp <= tolerance*norm_p)  
 #      if converged:  
 #          print "convergence after %s steps."%iter  
 #      else:  
 #          raise ArithmeticError("no convergence after %s steps."%iter)  
 #  
 #      return u,p  
             
 def MaskFromBoundaryTag(function_space,*tags):  
    """  
    create a mask on the given function space which one for samples  
    that touch the boundary tagged by tags.  
   
    usage: m=MaskFromBoundaryTag(Solution(domain),"left", "right")  
   
    @param function_space: a given function space  
    @type function_space: L{escript.FunctionSpace}  
    @param tags: boundray tags  
    @type tags: C{str}  
    @return: a mask which marks samples used by C{function_space} that are touching the  
             boundary tagged by any of the given tags.  
    @rtype: L{escript.Data} of rank 0  
    """  
    pde=linearPDEs.LinearPDE(function_space.getDomain(),numEquations=1, numSolutions=1)  
    d=escript.Scalar(0.,escript.FunctionOnBoundary(function_space.getDomain()))  
    for t in tags: d.setTaggedValue(t,1.)  
    pde.setValue(y=d)  
    out=util.whereNonZero(pde.getRightHandSide())  
    if out.getFunctionSpace() == function_space:  
       return out  
    else:  
       return util.whereNonZero(util.interpolate(out,function_space))  
   

Legend:
Removed from v.1312  
changed lines
  Added in v.1878

  ViewVC Help
Powered by ViewVC 1.1.26