/[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 1106 by gross, Thu Apr 19 01:54:57 2007 UTC revision 1956 by gross, Mon Nov 3 05:08:42 2008 UTC
# Line 1  Line 1 
1  # $Id$  
2    ########################################################
3    #
4    # Copyright (c) 2003-2008 by University of Queensland
5    # Earth Systems Science Computational Center (ESSCC)
6    # http://www.uq.edu.au/esscc
7    #
8    # Primary Business: Queensland, Australia
9    # Licensed under the Open Software License version 3.0
10    # http://www.opensource.org/licenses/osl-3.0.php
11    #
12    ########################################################
13    
14    __copyright__="""Copyright (c) 2003-2008 by University of Queensland
15    Earth Systems Science Computational Center (ESSCC)
16    http://www.uq.edu.au/esscc
17    Primary Business: Queensland, Australia"""
18    __license__="""Licensed under the Open Software License version 3.0
19    http://www.opensource.org/licenses/osl-3.0.php"""
20    __url__="http://www.uq.edu.au/esscc/escript-finley"
21    
22  """  """
23  Provides some tools related to PDEs.  Provides some tools related to PDEs.
# Line 18  Currently includes: Line 37  Currently includes:
37  """  """
38    
39  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
 __copyright__="""  Copyright (c) 2006 by ACcESS MNRF  
                     http://www.access.edu.au  
                 Primary Business: Queensland, Australia"""  
 __license__="""Licensed under the Open Software License version 3.0  
              http://www.opensource.org/licenses/osl-3.0.php"""  
 __url__="http://www.iservo.edu.au/esys"  
 __version__="$Revision$"  
 __date__="$Date$"  
40    
41    
42  import escript  import escript
43  import linearPDEs  import linearPDEs
44  import numarray  import numarray
45  import util  import util
46    import math
47    
48    ##### Added by Artak
49    # from Numeric import zeros,Int,Float64
50    ###################################
51    
52    
53  class TimeIntegrationManager:  class TimeIntegrationManager:
54    """    """
# Line 132  class Projector: Line 149  class Projector:
149      self.__pde.setValue(D = 1.)      self.__pde.setValue(D = 1.)
150      return      return
151    
   def __del__(self):  
     return  
   
152    def __call__(self, input_data):    def __call__(self, input_data):
153      """      """
154      Projects input_data onto a continuous function      Projects input_data onto a continuous function
# Line 142  class Projector: Line 156  class Projector:
156      @param input_data: The input_data to be projected.      @param input_data: The input_data to be projected.
157      """      """
158      out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())      out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())
159        self.__pde.setValue(Y = escript.Data(), Y_reduced = escript.Data())
160      if input_data.getRank()==0:      if input_data.getRank()==0:
161          self.__pde.setValue(Y = input_data)          self.__pde.setValue(Y = input_data)
162          out=self.__pde.getSolution()          out=self.__pde.getSolution()
# Line 395  class Locator: Line 410  class Locator:
410          else:          else:
411             return data             return data
412    
413    class SolverSchemeException(Exception):
414       """
415       exceptions thrown by solvers
416       """
417       pass
418    
419    class IndefinitePreconditioner(SolverSchemeException):
420       """
421       the preconditioner is not positive definite.
422       """
423       pass
424    class MaxIterReached(SolverSchemeException):
425       """
426       maxium number of iteration steps is reached.
427       """
428       pass
429    class IterationBreakDown(SolverSchemeException):
430       """
431       iteration scheme econouters an incurable breakdown.
432       """
433       pass
434    class NegativeNorm(SolverSchemeException):
435       """
436       a norm calculation returns a negative norm.
437       """
438       pass
439    
440    class IterationHistory(object):
441       """
442       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          @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          
465           @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           """
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    def PCG(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
497       """
498       Solver for
499    
500       M{Ax=b}
501    
502       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 iteration is terminated if the C{stoppingcriterium} function return C{True}.
506    
507       For details on the preconditioned conjugate gradient method see the book:
508    
509       Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
510       T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
511       C. Romine, and H. van der Vorst.
512    
513       @param b: the right hand side of the liner system. C{b} is altered.
514       @type b: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
515       @param Aprod: returns the value Ax
516       @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
518       @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>}
521       @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 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 stoppingcriterium: function that returns C{True} or C{False}
524       @param x: an initial guess for the solution. If no C{x} is given 0*b is used.
525       @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
526       @param iter_max: maximum number of iteration steps.
527       @type iter_max: C{int}
528       @return: the solution approximation and the corresponding residual
529       @rtype: C{tuple}
530       @warning: C{b} and C{x} are altered.
531       """
532       iter=0
533       if x==None:
534          x=0*b
535       else:
536          b += (-1)*Aprod(x)
537       r=b
538       rhat=Msolve(r)
539       d = rhat
540       rhat_dot_r = bilinearform(rhat, r)
541       if rhat_dot_r<0: raise NegativeNorm,"negative norm."
542    
543       while not stoppingcriterium(math.sqrt(rhat_dot_r),r,x):
544           iter+=1
545           if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
546    
547           q=Aprod(d)
548           alpha = rhat_dot_r / bilinearform(d, q)
549           x += alpha * d
550           r += (-alpha) * q
551    
552           rhat=Msolve(r)
553           rhat_dot_r_new = bilinearform(rhat, r)
554           beta = rhat_dot_r_new / rhat_dot_r
555           rhat+=beta * d
556           d=rhat
557    
558           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           try:  
1261               l=len(other)
1262               if l!=len(self):
1263                   raise ValueError,"length of of arguments don't match."
1264               for i in range(l): out.append(self[i]*other[i])
1265           except TypeError:
1266           for i in range(len(self)): out.append(self[i]*other)
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           try:  
1280               l=len(other)
1281               if l!=len(self):
1282                   raise ValueError,"length of of arguments don't match."
1283               for i in range(l): out.append(other[i]*self[i])
1284           except TypeError:
1285           for i in range(len(self)): out.append(other*self[i])
1286           return ArithmeticTuple(*tuple(out))
1287    
1288       def __div__(self,other):
1289           """
1290           dividing from the right
1291    
1292           @param other: scaling factor
1293           @type other: C{float}
1294           @return: itemwise self/other
1295           @rtype: L{ArithmeticTuple}
1296           """
1297           return self*(1/other)
1298    
1299       def __rdiv__(self,other):
1300           """
1301           dividing from the left
1302    
1303           @param other: scaling factor
1304           @type other: C{float}
1305           @return: itemwise other/self
1306           @rtype: L{ArithmeticTuple}
1307           """
1308           out=[]
1309           try:  
1310               l=len(other)
1311               if l!=len(self):
1312                   raise ValueError,"length of of arguments don't match."
1313               for i in range(l): out.append(other[i]/self[i])
1314           except TypeError:
1315           for i in range(len(self)): out.append(other/self[i])
1316           return ArithmeticTuple(*tuple(out))
1317      
1318       def __iadd__(self,other):
1319           """
1320           in-place add of other to self
1321    
1322           @param other: increment
1323           @type other: C{ArithmeticTuple}
1324           """
1325           if len(self) != len(other):
1326               raise ValueError,"tuple length must match."
1327           for i in range(len(self)):
1328               self.__items[i]+=other[i]
1329           return self
1330    
1331       def __add__(self,other):
1332           """
1333           add other to self
1334    
1335           @param other: increment
1336           @type other: C{ArithmeticTuple}
1337           """
1338           out=[]
1339           try:  
1340               l=len(other)
1341               if l!=len(self):
1342                   raise ValueError,"length of of arguments don't match."
1343               for i in range(l): out.append(self[i]+other[i])
1344           except TypeError:
1345           for i in range(len(self)): out.append(self[i]+other)
1346           return ArithmeticTuple(*tuple(out))
1347    
1348       def __sub__(self,other):
1349           """
1350           subtract other from self
1351    
1352           @param other: increment
1353           @type other: C{ArithmeticTuple}
1354           """
1355           out=[]
1356           try:  
1357               l=len(other)
1358               if l!=len(self):
1359                   raise ValueError,"length of of arguments don't match."
1360               for i in range(l): out.append(self[i]-other[i])
1361           except TypeError:
1362           for i in range(len(self)): out.append(self[i]-other)
1363           return ArithmeticTuple(*tuple(out))
1364      
1365       def __isub__(self,other):
1366           """
1367           subtract other from self
1368    
1369           @param other: increment
1370           @type other: C{ArithmeticTuple}
1371           """
1372           if len(self) != len(other):
1373               raise ValueError,"tuple length must match."
1374           for i in range(len(self)):
1375               self.__items[i]-=other[i]
1376           return self
1377    
1378       def __neg__(self):
1379           """
1380           negate
1381    
1382           """
1383           out=[]
1384           for i in range(len(self)):
1385               out.append(-self[i])
1386           return ArithmeticTuple(*tuple(out))
1387    
1388    
1389    class HomogeneousSaddlePointProblem(object):
1390          """
1391          This provides a framwork for solving linear homogeneous saddle point problem of the form
1392    
1393                 Av+B^*p=f
1394                 Bv    =0
1395    
1396          for the unknowns v and p and given operators A and B and given right hand side f.
1397          B^* is the adjoint operator of B is the given inner product.
1398    
1399          """
1400          def __init__(self,**kwargs):
1401            self.setTolerance()
1402            self.setToleranceReductionFactor()
1403    
1404          def initialize(self):
1405            """
1406            initialize the problem (overwrite)
1407            """
1408            pass
1409          def B(self,v):
1410             """
1411             returns Bv (overwrite)
1412             @rtype: equal to the type of p
1413    
1414             @note: boundary conditions on p should be zero!
1415             """
1416             pass
1417    
1418          def inner(self,p0,p1):
1419             """
1420             returns inner product of two element p0 and p1  (overwrite)
1421            
1422             @type p0: equal to the type of p
1423             @type p1: equal to the type of p
1424             @rtype: C{float}
1425    
1426             @rtype: equal to the type of p
1427             """
1428             pass
1429    
1430          def solve_A(self,u,p):
1431             """
1432             solves Av=f-Au-B^*p with accuracy self.getReducedTolerance() (overwrite)
1433    
1434             @rtype: equal to the type of v
1435             @note: boundary conditions on v should be zero!
1436             """
1437             pass
1438    
1439          def solve_prec(self,p):
1440             """
1441             provides a preconditioner for BA^{-1}B^* with accuracy self.getReducedTolerance() (overwrite)
1442    
1443             @rtype: equal to the type of p
1444             """
1445             pass
1446    
1447          def stoppingcriterium(self,Bv,v,p):
1448             """
1449             returns a True if iteration is terminated. (overwrite)
1450    
1451             @rtype: C{bool}
1452             """
1453             pass
1454                
1455          def __inner(self,p,r):
1456             return self.inner(p,r[1])
1457    
1458          def __inner_p(self,p1,p2):
1459             return self.inner(p1,p2)
1460          
1461          def __inner_a(self,a1,a2):
1462             return self.inner_a(a1,a2)
1463    
1464          def __inner_a1(self,a1,a2):
1465             return self.inner(a1[1],a2[1])
1466    
1467          def __stoppingcriterium(self,norm_r,r,p):
1468              return self.stoppingcriterium(r[1],r[0],p)
1469    
1470          def __stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None):
1471              return self.stoppingcriterium2(norm_r,norm_b,solver,TOL)
1472    
1473          def setTolerance(self,tolerance=1.e-8):
1474                  self.__tol=tolerance
1475          def getTolerance(self):
1476                  return self.__tol
1477          def setToleranceReductionFactor(self,reduction=0.01):
1478                  self.__reduction=reduction
1479          def getSubProblemTolerance(self):
1480                  return self.__reduction*self.getTolerance()
1481    
1482          def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='PCG',iter_restart=20):
1483                  """
1484                  solves the saddle point problem using initial guesses v and p.
1485    
1486                  @param max_iter: maximum number of iteration steps.
1487                  """
1488                  self.verbose=verbose
1489                  self.show_details=show_details and self.verbose
1490    
1491                  # assume p is known: then v=A^-1(f-B^*p)
1492                  # which leads to BA^-1B^*p = BA^-1f  
1493    
1494              # Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask)      
1495              self.__z=v+self.solve_A(v,p*0)
1496                  Bz=self.B(self.__z)
1497                  #
1498              #   solve BA^-1B^*p = Bz
1499                  #
1500                  #
1501                  #
1502                  self.iter=0
1503              if solver=='GMRES':      
1504                    if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter
1505                    p=GMRES(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.,iter_restart=iter_restart)
1506                    # solve Au=f-B^*p
1507                    #       A(u-v)=f-B^*p-Av
1508                    #       u=v+(u-v)
1509            u=v+self.solve_A(v,p)
1510    
1511              if solver=='TFQMR':      
1512                    if self.verbose: print "enter TFQMR method (iter_max=%s)"%max_iter
1513                    p=TFQMR(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)
1514                    # solve Au=f-B^*p
1515                    #       A(u-v)=f-B^*p-Av
1516                    #       u=v+(u-v)
1517            u=v+self.solve_A(v,p)
1518    
1519              if solver=='MINRES':      
1520                    if self.verbose: print "enter MINRES method (iter_max=%s)"%max_iter
1521                    p=MINRES(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)
1522                    # solve Au=f-B^*p
1523                    #       A(u-v)=f-B^*p-Av
1524                    #       u=v+(u-v)
1525            u=v+self.solve_A(v,p)
1526                  
1527              if solver=='GMRESC':      
1528                    if self.verbose: print "enter GMRES coupled method (iter_max=%s)"%max_iter
1529                    p0=self.solve_prec1(Bz)
1530                #alfa=(1/self.vol)*util.integrate(util.interpolate(p,escript.Function(self.domain)))
1531                    #p-=alfa
1532                    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)
1533                    #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())
1534    
1535                    # solve Au=f-B^*p
1536                    #       A(u-v)=f-B^*p-Av
1537                    #       u=v+(u-v)
1538                p=x[1]
1539            u=v+self.solve_A(v,p)      
1540            #p=x[1]
1541            #u=x[0]
1542    
1543                  if solver=='PCG':
1544                    #   note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv
1545                    #
1546                    #   with                    Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask)
1547                    #                           A(v-z)= f -Az - B^*p (v-z=0 on fixed_u_mask)
1548                    if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter
1549                    p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p)
1550                u=r[0]  
1551                    # print "DIFF=",util.integrate(p)
1552    
1553                  # print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)
1554    
1555              return u,p
1556    
1557          def __Msolve(self,r):
1558              return self.solve_prec(r[1])
1559    
1560          def __Msolve2(self,r):
1561              return self.solve_prec(r*1)
1562    
1563          def __Mempty(self,r):
1564              return r
1565    
1566    
1567          def __Aprod(self,p):
1568              # return BA^-1B*p
1569              #solve Av =B^*p as Av =f-Az-B^*(-p)
1570              v=self.solve_A(self.__z,-p)
1571              return ArithmeticTuple(v, self.B(v))
1572    
1573          def __Anext(self,x):
1574              # return next v,p
1575              #solve Av =-B^*p as Av =f-Az-B^*p
1576    
1577          pc=x[1]
1578              v=self.solve_A(self.__z,-pc)
1579          p=self.solve_prec1(self.B(v))
1580    
1581              return ArithmeticTuple(v,p)
1582    
1583    
1584          def __Aprod2(self,p):
1585              # return BA^-1B*p
1586              #solve Av =B^*p as Av =f-Az-B^*(-p)
1587          v=self.solve_A(self.__z,-p)
1588              return self.B(v)
1589    
1590          def __Aprod_Newton(self,p):
1591              # return BA^-1B*p - Bz
1592              #solve Av =-B^*p as Av =f-Az-B^*p
1593          v=self.solve_A(self.__z,-p)
1594              return self.B(v-self.__z)
1595    
1596          def __Aprod_Newton2(self,x):
1597              # return BA^-1B*p - Bz
1598              #solve Av =-B^*p as Av =f-Az-B^*p
1599              pc=x[1]
1600          v=self.solve_A(self.__z,-pc)
1601              p=self.solve_prec1(self.B(v-self.__z))
1602              return ArithmeticTuple(v,p)
1603    
1604    
1605    def MaskFromBoundaryTag(domain,*tags):
1606       """
1607       creates a mask on the Solution(domain) function space which one for samples
1608       that touch the boundary tagged by tags.
1609    
1610       usage: m=MaskFromBoundaryTag(domain,"left", "right")
1611    
1612       @param domain: a given domain
1613       @type domain: L{escript.Domain}
1614       @param tags: boundray tags
1615       @type tags: C{str}
1616       @return: a mask which marks samples that are touching the boundary tagged by any of the given tags.
1617       @rtype: L{escript.Data} of rank 0
1618       """
1619       pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1620       d=escript.Scalar(0.,escript.FunctionOnBoundary(domain))
1621       for t in tags: d.setTaggedValue(t,1.)
1622       pde.setValue(y=d)
1623       return util.whereNonZero(pde.getRightHandSide())
1624    #============================================================================================================================================
1625  class SaddlePointProblem(object):  class SaddlePointProblem(object):
1626     """     """
1627     This implements a solver for a saddlepoint problem     This implements a solver for a saddlepoint problem
# Line 420  class SaddlePointProblem(object): Line 1647  class SaddlePointProblem(object):
1647         @type verbose: C{bool}         @type verbose: C{bool}
1648         @note: this method may be overwritten by a particular saddle point problem         @note: this method may be overwritten by a particular saddle point problem
1649         """         """
1650           print "SaddlePointProblem should not be used anymore!"
1651           if not isinstance(verbose,bool):
1652                raise TypeError("verbose needs to be of type bool.")
1653         self.__verbose=verbose         self.__verbose=verbose
1654         self.relaxation=1.         self.relaxation=1.
1655           DeprecationWarning("SaddlePointProblem should not be used anymore.")
1656    
1657     def trace(self,text):     def trace(self,text):
1658         """         """
# Line 580  class SaddlePointProblem(object): Line 1811  class SaddlePointProblem(object):
1811              norm_u_new = util.Lsup(u_new)              norm_u_new = util.Lsup(u_new)
1812              p_new=p+self.relaxation*g_new              p_new=p+self.relaxation*g_new
1813              norm_p_new = util.sqrt(self.inner(p_new,p_new))              norm_p_new = util.sqrt(self.inner(p_new,p_new))
1814              self.trace("%s th step: f/u = %s, g/p = %s, relaxation = %s."%(self.iter,norm_f_new/norm_u_new, norm_g_new/norm_p_new, self.relaxation))              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))
1815    
1816              if self.iter>1:              if self.iter>1:
1817                 dg2=g_new-g                 dg2=g_new-g
# Line 595  class SaddlePointProblem(object): Line 1826  class SaddlePointProblem(object):
1826              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
1827          self.trace("convergence after %s steps."%self.iter)          self.trace("convergence after %s steps."%self.iter)
1828          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  
             
 # vim: expandtab shiftwidth=4:  

Legend:
Removed from v.1106  
changed lines
  Added in v.1956

  ViewVC Help
Powered by ViewVC 1.1.26