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

Legend:
Removed from v.877  
changed lines
  Added in v.1737

  ViewVC Help
Powered by ViewVC 1.1.26