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

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

  ViewVC Help
Powered by ViewVC 1.1.26