/[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 1484 by artak, Wed Apr 9 03:25:53 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):
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           self.history.append(norm_r)
494           if self.verbose: print "iter: %s:  norm(r) = %e"%(len(self.history)-1, self.history[-1])
495           return self.history[-1]<=self.tolerance * norm_b
496    
497    def PCG(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
498       """
499       Solver for
500    
501       M{Ax=b}
502    
503       with a symmetric and positive definite operator A (more details required!).
504       It uses the conjugate gradient method with preconditioner M providing an approximation of A.
505    
506       The iteration is terminated if the C{stoppingcriterium} function return C{True}.
507    
508       For details on the preconditioned conjugate gradient method see the book:
509    
510       Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
511       T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
512       C. Romine, and H. van der Vorst.
513    
514       @param b: the right hand side of the liner system. C{b} is altered.
515       @type b: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
516       @param Aprod: returns the value Ax
517       @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}.
518       @param Msolve: solves Mx=r
519       @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
520    type like argument C{x}.
521       @param bilinearform: inner product C{<x,r>}
522       @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}.
523       @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.
524       @type stoppingcriterium: function that returns C{True} or C{False}
525       @param x: an initial guess for the solution. If no C{x} is given 0*b is used.
526       @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
527       @param iter_max: maximum number of iteration steps.
528       @type iter_max: C{int}
529       @return: the solution approximation and the corresponding residual
530       @rtype: C{tuple}
531       @warning: C{b} and C{x} are altered.
532       """
533       iter=0
534       if x==None:
535          x=0*b
536       else:
537          b += (-1)*Aprod(x)
538       r=b
539       rhat=Msolve(r)
540       d = rhat
541       rhat_dot_r = bilinearform(rhat, r)
542       if rhat_dot_r<0: raise NegativeNorm,"negative norm."
543    
544       while not stoppingcriterium(math.sqrt(rhat_dot_r),r,x):
545           iter+=1
546           if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
547    
548           q=Aprod(d)
549           alpha = rhat_dot_r / bilinearform(d, q)
550           x += alpha * d
551           r += (-alpha) * q
552    
553           rhat=Msolve(r)
554           rhat_dot_r_new = bilinearform(rhat, r)
555           beta = rhat_dot_r_new / rhat_dot_r
556           rhat+=beta * d
557           d=rhat
558    
559           rhat_dot_r = rhat_dot_r_new
560           if rhat_dot_r<0: raise NegativeNorm,"negative norm."
561    
562       return x,r
563    
564    
565    ############################
566    # Added by Artak
567    #################################3
568    
569    #Apply a sequence of k Givens rotations, used within gmres codes
570    # vrot=givapp(c, s, vin, k)
571    def givapp(c,s,vin):
572        vrot=vin # warning: vin is altered!!!!
573        if isinstance(c,float):
574            vrot=[c*vrot[0]-s*vrot[1],s*vrot[0]+c*vrot[1]]
575        else:
576            for i in range(len(c)):
577                w1=c[i]*vrot[i]-s[i]*vrot[i+1]
578            w2=s[i]*vrot[i]+c[i]*vrot[i+1]
579                vrot[i:i+2]=w1,w2
580        return vrot
581    
582    def GMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=10):
583       m=iter_restart
584       iter=0
585       while True:
586          if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
587          x,stopped=GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=x, iter_max=iter_max-iter, iter_restart=m)
588          iter+=iter_restart
589          if stopped: break
590       return x
591    
592    def GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=10):
593       iter=0
594       r=Msolve(b)
595       r_dot_r = bilinearform(r, r)
596       if r_dot_r<0: raise NegativeNorm,"negative norm."
597       norm_b=math.sqrt(r_dot_r)
598    
599       if x==None:
600          x=0*b
601       else:
602          r=Msolve(b-Aprod(x))
603          r_dot_r = bilinearform(r, r)
604          if r_dot_r<0: raise NegativeNorm,"negative norm."
605      
606       h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)
607       c=numarray.zeros(iter_restart,numarray.Float64)
608       s=numarray.zeros(iter_restart,numarray.Float64)
609       g=numarray.zeros(iter_restart,numarray.Float64)
610       v=[]
611    
612       rho=math.sqrt(r_dot_r)
613       v.append(r/rho)
614       g[0]=rho
615    
616       while not (stoppingcriterium(rho,norm_b) or iter==iter_restart-1):
617    
618        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
619    
620        
621        p=Msolve(Aprod(v[iter]))
622    
623        v.append(p)
624    
625        v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
626    
627    # Modified Gram-Schmidt
628        for j in range(iter+1):
629          h[j][iter]=bilinearform(v[j],v[iter+1])  
630          v[iter+1]+=(-1.)*h[j][iter]*v[j]
631          
632        h[iter+1][iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
633        v_norm2=h[iter+1][iter]
634    
635    
636    # Reorthogonalize if needed
637        if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
638         for j in range(iter+1):
639            hr=bilinearform(v[j],v[iter+1])
640                h[j][iter]=h[j][iter]+hr #vhat
641                v[iter+1] +=(-1.)*hr*v[j]
642    
643         v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
644         h[iter+1][iter]=v_norm2
645    
646    #   watch out for happy breakdown
647            if v_norm2 != 0:
648             v[iter+1]=v[iter+1]/h[iter+1][iter]
649    
650    #   Form and store the information for the new Givens rotation
651        if iter > 0 :
652            hhat=[]
653            for i in range(iter+1) : hhat.append(h[i][iter])
654            hhat=givapp(c[0:iter],s[0:iter],hhat);
655                for i in range(iter+1) : h[i][iter]=hhat[i]
656    
657        mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])
658        if mu!=0 :
659            c[iter]=h[iter][iter]/mu
660            s[iter]=-h[iter+1][iter]/mu
661            h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter]
662            h[iter+1][iter]=0.0
663            g[iter:iter+2]=givapp(c[iter],s[iter],g[iter:iter+2])
664    
665    # Update the residual norm
666            rho=abs(g[iter+1])
667        iter+=1
668    
669    # At this point either iter > iter_max or rho < tol.
670    # It's time to compute x and leave.        
671    
672       if iter > 0 :
673         y=numarray.zeros(iter,numarray.Float64)    
674         y[iter-1] = g[iter-1] / h[iter-1][iter-1]
675         if iter > 1 :  
676            i=iter-2  
677            while i>=0 :
678              y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]
679              i=i-1
680         xhat=v[iter-1]*y[iter-1]
681         for i in range(iter-1):
682        xhat += v[i]*y[i]
683       else : xhat=v[0]
684        
685       x += xhat
686       if iter!=iter_restart-1:
687          stopped=True
688       else:
689          stopped=False
690    
691       return x,stopped
692    
693    def MINRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
694    
695        #
696        #  minres solves the system of linear equations Ax = b
697        #  where A is a symmetric matrix (possibly indefinite or singular)
698        #  and b is a given vector.
699        #  
700        #  "A" may be a dense or sparse matrix (preferably sparse!)
701        #  or the name of a function such that
702        #               y = A(x)
703        #  returns the product y = Ax for any given vector x.
704        #
705        #  "M" defines a positive-definite preconditioner M = C C'.
706        #  "M" may be a dense or sparse matrix (preferably sparse!)
707        #  or the name of a function such that
708        #  solves the system My = x for any given vector x.
709        #
710        #
711        
712        #------------------------------------------------------------------
713        # Set up y and v for the first Lanczos vector v1.
714        # y  =  beta1 P' v1,  where  P = C**(-1).
715        # v is really P' v1.
716        #------------------------------------------------------------------
717        if x==None:
718          x=0*b
719        else:
720          b += (-1)*Aprod(x)
721    
722        r1    = b
723        y = Msolve(b)
724        beta1 = bilinearform(b,y)
725    
726        if beta1< 0: raise NegativeNorm,"negative norm."
727    
728        #  If b = 0 exactly, stop with x = 0.
729        if beta1==0: return x*0.
730    
731        if beta1> 0:
732          beta1  = math.sqrt(beta1)       # Normalize y to get v1 later.
733    
734        #------------------------------------------------------------------
735        # Initialize quantities.
736        # ------------------------------------------------------------------
737        iter   = 0
738        Anorm = 0
739        ynorm = 0
740        oldb   = 0
741        beta   = beta1
742        dbar   = 0
743        epsln  = 0
744        phibar = beta1
745        rhs1   = beta1
746        rhs2   = 0
747        rnorm  = phibar
748        tnorm2 = 0
749        ynorm2 = 0
750        cs     = -1
751        sn     = 0
752        w      = b*0.
753        w2     = b*0.
754        r2     = r1
755        eps    = 0.0001
756    
757        #---------------------------------------------------------------------
758        # Main iteration loop.
759        # --------------------------------------------------------------------
760        while not stoppingcriterium(rnorm,Anorm*ynorm):    #  ||r|| / (||A|| ||x||)
761    
762        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
763            iter    = iter  +  1
764    
765            #-----------------------------------------------------------------
766            # Obtain quantities for the next Lanczos vector vk+1, k = 1, 2,...
767            # The general iteration is similar to the case k = 1 with v0 = 0:
768            #
769            #   p1      = Operator * v1  -  beta1 * v0,
770            #   alpha1  = v1'p1,
771            #   q2      = p2  -  alpha1 * v1,
772            #   beta2^2 = q2'q2,
773            #   v2      = (1/beta2) q2.
774            #
775            # Again, y = betak P vk,  where  P = C**(-1).
776            #-----------------------------------------------------------------
777            s = 1/beta                 # Normalize previous vector (in y).
778            v = s*y                    # v = vk if P = I
779        
780            y      = Aprod(v)
781        
782            if iter >= 2:
783              y = y - (beta/oldb)*r1
784    
785            alfa   = bilinearform(v,y)              # alphak
786            y      = (- alfa/beta)*r2 + y
787            r1     = r2
788            r2     = y
789            y = Msolve(r2)
790            oldb   = beta                           # oldb = betak
791            beta   = bilinearform(r2,y)             # beta = betak+1^2
792            if beta < 0: raise NegativeNorm,"negative norm."
793    
794            beta   = math.sqrt( beta )
795            tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta
796            
797            if iter==1:                 # Initialize a few things.
798              gmax   = abs( alfa )      # alpha1
799              gmin   = gmax             # alpha1
800    
801            # Apply previous rotation Qk-1 to get
802            #   [deltak epslnk+1] = [cs  sn][dbark    0   ]
803            #   [gbar k dbar k+1]   [sn -cs][alfak betak+1].
804        
805            oldeps = epsln
806            delta  = cs * dbar  +  sn * alfa  # delta1 = 0         deltak
807            gbar   = sn * dbar  -  cs * alfa  # gbar 1 = alfa1     gbar k
808            epsln  =               sn * beta  # epsln2 = 0         epslnk+1
809            dbar   =            -  cs * beta  # dbar 2 = beta2     dbar k+1
810    
811            # Compute the next plane rotation Qk
812    
813            gamma  = math.sqrt(gbar*gbar+beta*beta)  # gammak
814            gamma  = max(gamma,eps)
815            cs     = gbar / gamma             # ck
816            sn     = beta / gamma             # sk
817            phi    = cs * phibar              # phik
818            phibar = sn * phibar              # phibark+1
819    
820            # Update  x.
821    
822            denom = 1/gamma
823            w1    = w2
824            w2    = w
825            w     = (v - oldeps*w1 - delta*w2) * denom
826            x     = x  +  phi*w
827    
828            # Go round again.
829    
830            gmax   = max(gmax,gamma)
831            gmin   = min(gmin,gamma)
832            z      = rhs1 / gamma
833            ynorm2 = z*z  +  ynorm2
834            rhs1   = rhs2 -  delta*z
835            rhs2   =      -  epsln*z
836    
837            # Estimate various norms and test for convergence.
838    
839            Anorm  = math.sqrt( tnorm2 )
840            ynorm  = math.sqrt( ynorm2 )
841    
842            rnorm  = phibar
843    
844        return x
845        
846    #############################################
847    
848    class ArithmeticTuple(object):
849       """
850       tuple supporting inplace update x+=y and scaling x=a*y where x,y is an ArithmeticTuple and a is a float.
851    
852       example of usage:
853    
854       from esys.escript import Data
855       from numarray import array
856       a=Data(...)
857       b=array([1.,4.])
858       x=ArithmeticTuple(a,b)
859       y=5.*x
860    
861       """
862       def __init__(self,*args):
863           """
864           initialize object with elements args.
865    
866           @param args: tuple of object that support implace add (x+=y) and scaling (x=a*y)
867           """
868           self.__items=list(args)
869    
870       def __len__(self):
871           """
872           number of items
873    
874           @return: number of items
875           @rtype: C{int}
876           """
877           return len(self.__items)
878    
879       def __getitem__(self,index):
880           """
881           get an item
882    
883           @param index: item to be returned
884           @type index: C{int}
885           @return: item with index C{index}
886           """
887           return self.__items.__getitem__(index)
888    
889       def __mul__(self,other):
890           """
891           scaling from the right
892    
893           @param other: scaling factor
894           @type other: C{float}
895           @return: itemwise self*other
896           @rtype: L{ArithmeticTuple}
897           """
898           out=[]
899           for i in range(len(self)):
900               out.append(self[i]*other)
901           return ArithmeticTuple(*tuple(out))
902    
903       def __rmul__(self,other):
904           """
905           scaling from the left
906    
907           @param other: scaling factor
908           @type other: C{float}
909           @return: itemwise other*self
910           @rtype: L{ArithmeticTuple}
911           """
912           out=[]
913           for i in range(len(self)):
914               out.append(other*self[i])
915           return ArithmeticTuple(*tuple(out))
916    
917    #########################
918    # Added by Artak
919    #########################
920       def __div__(self,other):
921           """
922           dividing from the right
923    
924           @param other: scaling factor
925           @type other: C{float}
926           @return: itemwise self/other
927           @rtype: L{ArithmeticTuple}
928           """
929           out=[]
930           for i in range(len(self)):
931               out.append(self[i]/other)
932           return ArithmeticTuple(*tuple(out))
933    
934       def __rdiv__(self,other):
935           """
936           dividing from the left
937    
938           @param other: scaling factor
939           @type other: C{float}
940           @return: itemwise other/self
941           @rtype: L{ArithmeticTuple}
942           """
943           out=[]
944           for i in range(len(self)):
945               out.append(other/self[i])
946           return ArithmeticTuple(*tuple(out))
947      
948    ##########################################33
949    
950       def __iadd__(self,other):
951           """
952           in-place add of other to self
953    
954           @param other: increment
955           @type other: C{ArithmeticTuple}
956           """
957           if len(self) != len(other):
958               raise ValueError,"tuple length must match."
959           for i in range(len(self)):
960               self.__items[i]+=other[i]
961           return self
962    
963    class HomogeneousSaddlePointProblem(object):
964          """
965          This provides a framwork for solving homogeneous saddle point problem of the form
966    
967                 Av+B^*p=f
968                 Bv    =0
969    
970          for the unknowns v and p and given operators A and B and given right hand side f.
971          B^* is the adjoint operator of B is the given inner product.
972    
973          """
974          def __init__(self,**kwargs):
975            self.setTolerance()
976            self.setToleranceReductionFactor()
977    
978          def initialize(self):
979            """
980            initialize the problem (overwrite)
981            """
982            pass
983          def B(self,v):
984             """
985             returns Bv (overwrite)
986             @rtype: equal to the type of p
987    
988             @note: boundary conditions on p should be zero!
989             """
990             pass
991    
992          def inner(self,p0,p1):
993             """
994             returns inner product of two element p0 and p1  (overwrite)
995            
996             @type p0: equal to the type of p
997             @type p1: equal to the type of p
998             @rtype: C{float}
999    
1000             @rtype: equal to the type of p
1001             """
1002             pass
1003    
1004          def solve_A(self,u,p):
1005             """
1006             solves Av=f-Au-B^*p with accuracy self.getReducedTolerance() (overwrite)
1007    
1008             @rtype: equal to the type of v
1009             @note: boundary conditions on v should be zero!
1010             """
1011             pass
1012    
1013          def solve_prec(self,p):
1014             """
1015             provides a preconditioner for BA^{-1}B^* with accuracy self.getReducedTolerance() (overwrite)
1016    
1017             @rtype: equal to the type of p
1018             """
1019             pass
1020    
1021          def stoppingcriterium(self,Bv,v,p):
1022             """
1023             returns a True if iteration is terminated. (overwrite)
1024    
1025             @rtype: C{bool}
1026             """
1027             pass
1028                
1029          def __inner(self,p,r):
1030             return self.inner(p,r[1])
1031    
1032          def __inner_p(self,p1,p2):
1033             return self.inner(p1,p2)
1034    
1035          def __stoppingcriterium(self,norm_r,r,p):
1036              return self.stoppingcriterium(r[1],r[0],p)
1037    
1038          def __stoppingcriterium_GMRES(self,norm_r,norm_b):
1039              return self.stoppingcriterium_GMRES(norm_r,norm_b)
1040    
1041          def __stoppingcriterium_MINRES(self,norm_r,norm_Ax):
1042              return self.stoppingcriterium_MINRES(norm_r,norm_Ax)
1043    
1044    
1045          def setTolerance(self,tolerance=1.e-8):
1046                  self.__tol=tolerance
1047          def getTolerance(self):
1048                  return self.__tol
1049          def setToleranceReductionFactor(self,reduction=0.01):
1050                  self.__reduction=reduction
1051          def getSubProblemTolerance(self):
1052                  return self.__reduction*self.getTolerance()
1053    
1054          def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='PCG'):
1055                  """
1056                  solves the saddle point problem using initial guesses v and p.
1057    
1058                  @param max_iter: maximum number of iteration steps.
1059                  """
1060                  self.verbose=verbose
1061                  self.show_details=show_details and self.verbose
1062    
1063                  # assume p is known: then v=A^-1(f-B^*p)
1064                  # which leads to BA^-1B^*p = BA^-1f  
1065    
1066              # Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask)      
1067    
1068              
1069              self.__z=v+self.solve_A(v,p*0)
1070    
1071                  Bz=self.B(self.__z)
1072                  #
1073              #   solve BA^-1B^*p = Bz
1074                  #
1075                  #   note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv
1076                  #
1077                  #   with                    Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask)
1078                  #                           A(v-z)=Az-B^*p-Az = f -Az - B^*p (v-z=0 on fixed_u_mask)
1079                  #
1080                  self.iter=0
1081              if solver=='GMRES':      
1082                    if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter
1083                    p=GMRES(Bz,self.__Aprod_GMRES,self.__Msolve_GMRES,self.__inner_p,self.__stoppingcriterium_GMRES,iter_max=max_iter, x=p*1.)
1084                    # solve Au=f-B^*p
1085                    #       A(u-v)=f-B^*p-Av
1086                    #       u=v+(u-v)
1087            u=v+self.solve_A(v,p)
1088    
1089              if solver=='MINRES':      
1090                    if self.verbose: print "enter MINRES method (iter_max=%s)"%max_iter
1091                    p=GMRES(Bz,self.__Aprod_GMRES,self.__Msolve_GMRES,self.__inner_p,self.__stoppingcriterium_MINRES,iter_max=max_iter, x=p*1.)
1092                    # solve Au=f-B^*p
1093                    #       A(u-v)=f-B^*p-Av
1094                    #       u=v+(u-v)
1095            u=v+self.solve_A(v,p)
1096                  else:
1097                    if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter
1098                    p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p)
1099                u=r[0]  
1100    
1101                  print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)
1102    
1103              return u,p
1104    
1105          def __Msolve(self,r):
1106              return self.solve_prec(r[1])
1107    
1108          def __Msolve_GMRES(self,r):
1109              return self.solve_prec(r)
1110    
1111    
1112          def __Aprod(self,p):
1113              # return BA^-1B*p
1114              #solve Av =-B^*p as Av =f-Az-B^*p
1115              v=self.solve_A(self.__z,-p)
1116              return ArithmeticTuple(v, self.B(v))
1117    
1118          def __Aprod_GMRES(self,p):
1119              # return BA^-1B*p
1120              #solve Av =-B^*p as Av =f-Az-B^*p
1121          v=self.solve_A(self.__z,-p)
1122              return self.B(v)
1123    
1124  class SaddlePointProblem(object):  class SaddlePointProblem(object):
1125     """     """
1126     This implements a solver for a saddlepoint problem     This implements a solver for a saddlepoint problem
# Line 363  class SaddlePointProblem(object): Line 1130  class SaddlePointProblem(object):
1130    
1131     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:
1132    
1133     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})}
1134     M{Q_g (p^{k+1}-p^{k}) =   g(u^{k+1})}     M{Q_g (p^{k+1}-p^{k}) =   g(u^{k+1})}
1135    
1136     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 1142  class SaddlePointProblem(object):
1142         """         """
1143         initializes the problem         initializes the problem
1144    
1145         @parm verbose: switches on the printing out some information         @param verbose: switches on the printing out some information
1146         @type verbose: C{bool}         @type verbose: C{bool}
1147         @note: this method may be overwritten by a particular saddle point problem         @note: this method may be overwritten by a particular saddle point problem
1148         """         """
1149           if not isinstance(verbose,bool):
1150                raise TypeError("verbose needs to be of type bool.")
1151         self.__verbose=verbose         self.__verbose=verbose
1152         self.relaxation=1.         self.relaxation=1.
1153    
# Line 386  class SaddlePointProblem(object): Line 1155  class SaddlePointProblem(object):
1155         """         """
1156         prints text if verbose has been set         prints text if verbose has been set
1157    
1158         @parm text: a text message         @param text: a text message
1159         @type text: C{str}         @type text: C{str}
1160         """         """
1161         if self.__verbose: print "%s: %s"%(str(self),text)         if self.__verbose: print "%s: %s"%(str(self),text)
# Line 539  class SaddlePointProblem(object): Line 1308  class SaddlePointProblem(object):
1308              norm_u_new = util.Lsup(u_new)              norm_u_new = util.Lsup(u_new)
1309              p_new=p+self.relaxation*g_new              p_new=p+self.relaxation*g_new
1310              norm_p_new = util.sqrt(self.inner(p_new,p_new))              norm_p_new = util.sqrt(self.inner(p_new,p_new))
1311              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))
1312    
1313              if self.iter>1:              if self.iter>1:
1314                 dg2=g_new-g                 dg2=g_new-g
# Line 551  class SaddlePointProblem(object): Line 1320  class SaddlePointProblem(object):
1320                 tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new                 tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new
1321                 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:
1322                     converged=True                     converged=True
                    break  
1323              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
1324          self.trace("convergence after %s steps."%self.iter)          self.trace("convergence after %s steps."%self.iter)
1325          return u,p          return u,p
# Line 582  class SaddlePointProblem(object): Line 1350  class SaddlePointProblem(object):
1350  #  #
1351  #      return u,p  #      return u,p
1352                        
1353  # vim: expandtab shiftwidth=4:  def MaskFromBoundaryTag(function_space,*tags):
1354       """
1355       create a mask on the given function space which one for samples
1356       that touch the boundary tagged by tags.
1357    
1358       usage: m=MaskFromBoundaryTag(Solution(domain),"left", "right")
1359    
1360       @param function_space: a given function space
1361       @type function_space: L{escript.FunctionSpace}
1362       @param tags: boundray tags
1363       @type tags: C{str}
1364       @return: a mask which marks samples used by C{function_space} that are touching the
1365                boundary tagged by any of the given tags.
1366       @rtype: L{escript.Data} of rank 0
1367       """
1368       pde=linearPDEs.LinearPDE(function_space.getDomain(),numEquations=1, numSolutions=1)
1369       d=escript.Scalar(0.,escript.FunctionOnBoundary(function_space.getDomain()))
1370       for t in tags: d.setTaggedValue(t,1.)
1371       pde.setValue(y=d)
1372       out=util.whereNonZero(pde.getRightHandSide())
1373       if out.getFunctionSpace() == function_space:
1374          return out
1375       else:
1376          return util.whereNonZero(util.interpolate(out,function_space))
1377    
1378    
1379    

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

  ViewVC Help
Powered by ViewVC 1.1.26