/[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 782 by bcumming, Tue Jul 18 00:47:47 2006 UTC revision 1482 by artak, Wed Apr 9 02:29:47 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 7  Currently includes: Line 21  Currently includes:
21      - Projector - to project a discontinuous      - Projector - to project a discontinuous
22      - Locator - to trace values in data objects at a certain location      - Locator - to trace values in data objects at a certain location
23      - TimeIntegrationManager - to handel extraplotion in time      - TimeIntegrationManager - to handel extraplotion in time
24        - SaddlePointProblem - solver for Saddle point problems using the inexact uszawa scheme
25    
26  @var __author__: name of author  @var __author__: name of author
27  @var __copyright__: copyrights  @var __copyright__: copyrights
# Line 31  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 131  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 141  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 296  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 315  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 341  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  # vim: expandtab shiftwidth=4:  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 other quantities.
736        # ------------------------------------------------------------------
737    #  Initialize                              
738    
739        iter   = 0
740        Anorm = 0
741        ynorm = 0
742    #    x=x*0
743    
744        oldb   = 0
745        beta   = beta1
746        dbar   = 0
747        epsln  = 0
748        phibar = beta1
749        rhs1   = beta1
750        rhs2   = 0
751        rnorm  = phibar
752        tnorm2 = 0
753        ynorm2 = 0
754        cs     = -1
755        sn     = 0
756        w      = b*0.
757        w2     = b*0.
758        r2     = r1
759        eps    = 0.0001
760    
761        #---------------------------------------------------------------------
762        # Main iteration loop.
763        # --------------------------------------------------------------------
764        while not stoppingcriterium(rnorm,Anorm*ynorm):    #  ||r|| / (||A|| ||x||)
765    
766        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
767            iter    = iter  +  1
768    
769            #-----------------------------------------------------------------
770            # Obtain quantities for the next Lanczos vector vk+1, k = 1, 2,...
771            # The general iteration is similar to the case k = 1 with v0 = 0:
772            #
773            #   p1      = Operator * v1  -  beta1 * v0,
774            #   alpha1  = v1'p1,
775            #   q2      = p2  -  alpha1 * v1,
776            #   beta2^2 = q2'q2,
777            #   v2      = (1/beta2) q2.
778            #
779            # Again, y = betak P vk,  where  P = C**(-1).
780            #-----------------------------------------------------------------
781            s = 1/beta                 # Normalize previous vector (in y).
782            v = s*y                    # v = vk if P = I
783        
784            y      = Aprod(v)
785        
786            if iter >= 2:
787              y = y - (beta/oldb)*r1
788    
789            alfa   = bilinearform(v,y)              # alphak
790            y      = (- alfa/beta)*r2 + y
791            r1     = r2
792            r2     = y
793            y = Msolve(r2)
794            oldb   = beta                           # oldb = betak
795            beta   = bilinearform(r2,y)             # beta = betak+1^2
796            if beta < 0: raise NegativeNorm,"negative norm."
797    
798            beta   = math.sqrt( beta )
799            tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta
800            
801            if iter==1:                 # Initialize a few things.
802              gmax   = abs( alfa )      # alpha1
803              gmin   = gmax             # alpha1
804    
805            # Apply previous rotation Qk-1 to get
806            #   [deltak epslnk+1] = [cs  sn][dbark    0   ]
807            #   [gbar k dbar k+1]   [sn -cs][alfak betak+1].
808        
809            oldeps = epsln
810            delta  = cs * dbar  +  sn * alfa  # delta1 = 0         deltak
811            gbar   = sn * dbar  -  cs * alfa  # gbar 1 = alfa1     gbar k
812            epsln  =               sn * beta  # epsln2 = 0         epslnk+1
813            dbar   =            -  cs * beta  # dbar 2 = beta2     dbar k+1
814    
815            # Compute the next plane rotation Qk
816    
817            gamma  = math.sqrt(gbar*gbar+beta*beta)  # gammak
818            gamma  = max(gamma,eps)
819            cs     = gbar / gamma             # ck
820            sn     = beta / gamma             # sk
821            phi    = cs * phibar              # phik
822            phibar = sn * phibar              # phibark+1
823    
824            # Update  x.
825    
826            denom = 1/gamma
827            w1    = w2
828            w2    = w
829            w     = (v - oldeps*w1 - delta*w2) * denom
830            x     = x  +  phi*w
831    
832            # Go round again.
833    
834            gmax   = max(gmax,gamma)
835            gmin   = min(gmin,gamma)
836            z      = rhs1 / gamma
837            ynorm2 = z*z  +  ynorm2
838            rhs1   = rhs2 -  delta*z
839            rhs2   =      -  epsln*z
840    
841            # Estimate various norms and test for convergence.
842    
843            Anorm  = math.sqrt( tnorm2 )
844            ynorm  = math.sqrt( ynorm2 )
845    
846            rnorm  = phibar
847    
848        # Return final answer.
849        print iter
850        return x
851        
852    #############################################
853    
854    class ArithmeticTuple(object):
855       """
856       tuple supporting inplace update x+=y and scaling x=a*y where x,y is an ArithmeticTuple and a is a float.
857    
858       example of usage:
859    
860       from esys.escript import Data
861       from numarray import array
862       a=Data(...)
863       b=array([1.,4.])
864       x=ArithmeticTuple(a,b)
865       y=5.*x
866    
867       """
868       def __init__(self,*args):
869           """
870           initialize object with elements args.
871    
872           @param args: tuple of object that support implace add (x+=y) and scaling (x=a*y)
873           """
874           self.__items=list(args)
875    
876       def __len__(self):
877           """
878           number of items
879    
880           @return: number of items
881           @rtype: C{int}
882           """
883           return len(self.__items)
884    
885       def __getitem__(self,index):
886           """
887           get an item
888    
889           @param index: item to be returned
890           @type index: C{int}
891           @return: item with index C{index}
892           """
893           return self.__items.__getitem__(index)
894    
895       def __mul__(self,other):
896           """
897           scaling from the right
898    
899           @param other: scaling factor
900           @type other: C{float}
901           @return: itemwise self*other
902           @rtype: L{ArithmeticTuple}
903           """
904           out=[]
905           for i in range(len(self)):
906               out.append(self[i]*other)
907           return ArithmeticTuple(*tuple(out))
908    
909       def __rmul__(self,other):
910           """
911           scaling from the left
912    
913           @param other: scaling factor
914           @type other: C{float}
915           @return: itemwise other*self
916           @rtype: L{ArithmeticTuple}
917           """
918           out=[]
919           for i in range(len(self)):
920               out.append(other*self[i])
921           return ArithmeticTuple(*tuple(out))
922    
923    #########################
924    # Added by Artak
925    #########################
926       def __div__(self,other):
927           """
928           dividing from the right
929    
930           @param other: scaling factor
931           @type other: C{float}
932           @return: itemwise self/other
933           @rtype: L{ArithmeticTuple}
934           """
935           out=[]
936           for i in range(len(self)):
937               out.append(self[i]/other)
938           return ArithmeticTuple(*tuple(out))
939    
940       def __rdiv__(self,other):
941           """
942           dividing from the left
943    
944           @param other: scaling factor
945           @type other: C{float}
946           @return: itemwise other/self
947           @rtype: L{ArithmeticTuple}
948           """
949           out=[]
950           for i in range(len(self)):
951               out.append(other/self[i])
952           return ArithmeticTuple(*tuple(out))
953      
954    ##########################################33
955    
956       def __iadd__(self,other):
957           """
958           in-place add of other to self
959    
960           @param other: increment
961           @type other: C{ArithmeticTuple}
962           """
963           if len(self) != len(other):
964               raise ValueError,"tuple length must match."
965           for i in range(len(self)):
966               self.__items[i]+=other[i]
967           return self
968    
969    class HomogeneousSaddlePointProblem(object):
970          """
971          This provides a framwork for solving homogeneous saddle point problem of the form
972    
973                 Av+B^*p=f
974                 Bv    =0
975    
976          for the unknowns v and p and given operators A and B and given right hand side f.
977          B^* is the adjoint operator of B is the given inner product.
978    
979          """
980          def __init__(self,**kwargs):
981            self.setTolerance()
982            self.setToleranceReductionFactor()
983    
984          def initialize(self):
985            """
986            initialize the problem (overwrite)
987            """
988            pass
989          def B(self,v):
990             """
991             returns Bv (overwrite)
992             @rtype: equal to the type of p
993    
994             @note: boundary conditions on p should be zero!
995             """
996             pass
997    
998          def inner(self,p0,p1):
999             """
1000             returns inner product of two element p0 and p1  (overwrite)
1001            
1002             @type p0: equal to the type of p
1003             @type p1: equal to the type of p
1004             @rtype: C{float}
1005    
1006             @rtype: equal to the type of p
1007             """
1008             pass
1009    
1010          def solve_A(self,u,p):
1011             """
1012             solves Av=f-Au-B^*p with accuracy self.getReducedTolerance() (overwrite)
1013    
1014             @rtype: equal to the type of v
1015             @note: boundary conditions on v should be zero!
1016             """
1017             pass
1018    
1019          def solve_prec(self,p):
1020             """
1021             provides a preconditioner for BA^{-1}B^* with accuracy self.getReducedTolerance() (overwrite)
1022    
1023             @rtype: equal to the type of p
1024             """
1025             pass
1026    
1027          def stoppingcriterium(self,Bv,v,p):
1028             """
1029             returns a True if iteration is terminated. (overwrite)
1030    
1031             @rtype: C{bool}
1032             """
1033             pass
1034                
1035          def __inner(self,p,r):
1036             return self.inner(p,r[1])
1037    
1038          def __inner_p(self,p1,p2):
1039             return self.inner(p1,p2)
1040    
1041          def __stoppingcriterium(self,norm_r,r,p):
1042              return self.stoppingcriterium(r[1],r[0],p)
1043    
1044          def __stoppingcriterium_GMRES(self,norm_r,norm_b):
1045              return self.stoppingcriterium_GMRES(norm_r,norm_b)
1046    
1047          def __stoppingcriterium_MINRES(self,norm_r,norm_Ax):
1048              return self.stoppingcriterium_MINRES(norm_r,norm_Ax)
1049    
1050    
1051          def setTolerance(self,tolerance=1.e-8):
1052                  self.__tol=tolerance
1053          def getTolerance(self):
1054                  return self.__tol
1055          def setToleranceReductionFactor(self,reduction=0.01):
1056                  self.__reduction=reduction
1057          def getSubProblemTolerance(self):
1058                  return self.__reduction*self.getTolerance()
1059    
1060          def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='PCG'):
1061                  """
1062                  solves the saddle point problem using initial guesses v and p.
1063    
1064                  @param max_iter: maximum number of iteration steps.
1065                  """
1066                  self.verbose=verbose
1067                  self.show_details=show_details and self.verbose
1068    
1069                  # assume p is known: then v=A^-1(f-B^*p)
1070                  # which leads to BA^-1B^*p = BA^-1f  
1071    
1072              # Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask)      
1073    
1074              
1075              self.__z=v+self.solve_A(v,p*0)
1076    
1077                  Bz=self.B(self.__z)
1078                  #
1079              #   solve BA^-1B^*p = Bz
1080                  #
1081                  #   note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv
1082                  #
1083                  #   with                    Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask)
1084                  #                           A(v-z)=Az-B^*p-Az = f -Az - B^*p (v-z=0 on fixed_u_mask)
1085                  #
1086                  self.iter=0
1087              if solver=='GMRES':      
1088                    if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter
1089                    p=GMRES(Bz,self.__Aprod_GMRES,self.__Msolve_GMRES,self.__inner_p,self.__stoppingcriterium_GMRES,iter_max=max_iter, x=p*1.)
1090                    # solve Au=f-B^*p
1091                    #       A(u-v)=f-B^*p-Av
1092                    #       u=v+(u-v)
1093            u=v+self.solve_A(v,p)
1094    
1095              if solver=='MINRES':      
1096                    if self.verbose: print "enter MINRES method (iter_max=%s)"%max_iter
1097                    p=GMRES(Bz,self.__Aprod_GMRES,self.__Msolve_GMRES,self.__inner_p,self.__stoppingcriterium_MINRES,iter_max=max_iter, x=p*1.)
1098                    # solve Au=f-B^*p
1099                    #       A(u-v)=f-B^*p-Av
1100                    #       u=v+(u-v)
1101            u=v+self.solve_A(v,p)
1102        
1103                  if solver=='PCG':
1104                    if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter
1105                    p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p)
1106                u=r[0]  
1107    
1108                  print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)
1109    
1110              return u,p
1111    
1112          def __Msolve(self,r):
1113              return self.solve_prec(r[1])
1114    
1115          def __Msolve_GMRES(self,r):
1116              return self.solve_prec(r)
1117    
1118    
1119          def __Aprod(self,p):
1120              # return BA^-1B*p
1121              #solve Av =-B^*p as Av =f-Az-B^*p
1122              v=self.solve_A(self.__z,-p)
1123              return ArithmeticTuple(v, self.B(v))
1124    
1125          def __Aprod_GMRES(self,p):
1126              # return BA^-1B*p
1127              #solve Av =-B^*p as Av =f-Az-B^*p
1128          v=self.solve_A(self.__z,-p)
1129              return self.B(v)
1130    
1131    class SaddlePointProblem(object):
1132       """
1133       This implements a solver for a saddlepoint problem
1134    
1135       M{f(u,p)=0}
1136       M{g(u)=0}
1137    
1138       for u and p. The problem is solved with an inexact Uszawa scheme for p:
1139    
1140       M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})}
1141       M{Q_g (p^{k+1}-p^{k}) =   g(u^{k+1})}
1142    
1143       where Q_f is an approximation of the Jacobiean A_f of f with respect to u  and Q_f is an approximation of
1144       A_g A_f^{-1} A_g with A_g is the jacobiean of g with respect to p. As a the construction of a 'proper'
1145       Q_g can be difficult, non-linear conjugate gradient method is applied to solve for p, so Q_g plays
1146       in fact the role of a preconditioner.
1147       """
1148       def __init__(self,verbose=False,*args):
1149           """
1150           initializes the problem
1151    
1152           @param verbose: switches on the printing out some information
1153           @type verbose: C{bool}
1154           @note: this method may be overwritten by a particular saddle point problem
1155           """
1156           if not isinstance(verbose,bool):
1157                raise TypeError("verbose needs to be of type bool.")
1158           self.__verbose=verbose
1159           self.relaxation=1.
1160    
1161       def trace(self,text):
1162           """
1163           prints text if verbose has been set
1164    
1165           @param text: a text message
1166           @type text: C{str}
1167           """
1168           if self.__verbose: print "%s: %s"%(str(self),text)
1169    
1170       def solve_f(self,u,p,tol=1.e-8):
1171           """
1172           solves
1173    
1174           A_f du = f(u,p)
1175    
1176           with tolerance C{tol} and return du. A_f is Jacobiean of f with respect to u.
1177    
1178           @param u: current approximation of u
1179           @type u: L{escript.Data}
1180           @param p: current approximation of p
1181           @type p: L{escript.Data}
1182           @param tol: tolerance expected for du
1183           @type tol: C{float}
1184           @return: increment du
1185           @rtype: L{escript.Data}
1186           @note: this method has to be overwritten by a particular saddle point problem
1187           """
1188           pass
1189    
1190       def solve_g(self,u,tol=1.e-8):
1191           """
1192           solves
1193    
1194           Q_g dp = g(u)
1195    
1196           with Q_g is a preconditioner for A_g A_f^{-1} A_g with  A_g is the jacobiean of g with respect to p.
1197    
1198           @param u: current approximation of u
1199           @type u: L{escript.Data}
1200           @param tol: tolerance expected for dp
1201           @type tol: C{float}
1202           @return: increment dp
1203           @rtype: L{escript.Data}
1204           @note: this method has to be overwritten by a particular saddle point problem
1205           """
1206           pass
1207    
1208       def inner(self,p0,p1):
1209           """
1210           inner product of p0 and p1 approximating p. Typically this returns integrate(p0*p1)
1211           @return: inner product of p0 and p1
1212           @rtype: C{float}
1213           """
1214           pass
1215    
1216       subiter_max=3
1217       def solve(self,u0,p0,tolerance=1.e-6,tolerance_u=None,iter_max=100,accepted_reduction=0.995,relaxation=None):
1218            """
1219            runs the solver
1220    
1221            @param u0: initial guess for C{u}
1222            @type u0: L{esys.escript.Data}
1223            @param p0: initial guess for C{p}
1224            @type p0: L{esys.escript.Data}
1225            @param tolerance: tolerance for relative error in C{u} and C{p}
1226            @type tolerance: positive C{float}
1227            @param tolerance_u: tolerance for relative error in C{u} if different from C{tolerance}
1228            @type tolerance_u: positive C{float}
1229            @param iter_max: maximum number of iteration steps.
1230            @type iter_max: C{int}
1231            @param accepted_reduction: if the norm  g cannot be reduced by C{accepted_reduction} backtracking to adapt the
1232                                       relaxation factor. If C{accepted_reduction=None} no backtracking is used.
1233            @type accepted_reduction: positive C{float} or C{None}
1234            @param relaxation: initial relaxation factor. If C{relaxation==None}, the last relaxation factor is used.
1235            @type relaxation: C{float} or C{None}
1236            """
1237            tol=1.e-2
1238            if tolerance_u==None: tolerance_u=tolerance
1239            if not relaxation==None: self.relaxation=relaxation
1240            if accepted_reduction ==None:
1241                  angle_limit=0.
1242            elif accepted_reduction>=1.:
1243                  angle_limit=0.
1244            else:
1245                  angle_limit=util.sqrt(1-accepted_reduction**2)
1246            self.iter=0
1247            u=u0
1248            p=p0
1249            #
1250            #   initialize things:
1251            #
1252            converged=False
1253            #
1254            #  start loop:
1255            #
1256            #  initial search direction is g
1257            #
1258            while not converged :
1259                if self.iter>iter_max:
1260                    raise ArithmeticError("no convergence after %s steps."%self.iter)
1261                f_new=self.solve_f(u,p,tol)
1262                norm_f_new = util.Lsup(f_new)
1263                u_new=u-f_new
1264                g_new=self.solve_g(u_new,tol)
1265                self.iter+=1
1266                norm_g_new = util.sqrt(self.inner(g_new,g_new))
1267                if norm_f_new==0. and norm_g_new==0.: return u, p
1268                if self.iter>1 and not accepted_reduction==None:
1269                   #
1270                   #   did we manage to reduce the norm of G? I
1271                   #   if not we start a backtracking procedure
1272                   #
1273                   # print "new/old norm = ",norm_g_new, norm_g, norm_g_new/norm_g
1274                   if norm_g_new > accepted_reduction * norm_g:
1275                      sub_iter=0
1276                      s=self.relaxation
1277                      d=g
1278                      g_last=g
1279                      self.trace("    start substepping: f = %s, g = %s, relaxation = %s."%(norm_f_new, norm_g_new, s))
1280                      while sub_iter < self.subiter_max and  norm_g_new > accepted_reduction * norm_g:
1281                         dg= g_new-g_last
1282                         norm_dg=abs(util.sqrt(self.inner(dg,dg))/self.relaxation)
1283                         rad=self.inner(g_new,dg)/self.relaxation
1284                         # print "   ",sub_iter,": rad, norm_dg:",abs(rad), norm_dg*norm_g_new * angle_limit
1285                         # print "   ",sub_iter,": rad, norm_dg:",rad, norm_dg, norm_g_new, norm_g
1286                         if abs(rad) < norm_dg*norm_g_new * angle_limit:
1287                             if sub_iter>0: self.trace("    no further improvements expected from backtracking.")
1288                             break
1289                         r=self.relaxation
1290                         self.relaxation= - rad/norm_dg**2
1291                         s+=self.relaxation
1292                         #####
1293                         # a=g_new+self.relaxation*dg/r
1294                         # print "predicted new norm = ",util.sqrt(self.inner(a,a)),util.sqrt(self.inner(g_new,g_new)), self.relaxation
1295                         #####
1296                         g_last=g_new
1297                         p+=self.relaxation*d
1298                         f_new=self.solve_f(u,p,tol)
1299                         u_new=u-f_new
1300                         g_new=self.solve_g(u_new,tol)
1301                         self.iter+=1
1302                         norm_f_new = util.Lsup(f_new)
1303                         norm_g_new = util.sqrt(self.inner(g_new,g_new))
1304                         # print "   ",sub_iter," new g norm",norm_g_new
1305                         self.trace("    %s th sub-step: f = %s, g = %s, relaxation = %s."%(sub_iter, norm_f_new, norm_g_new, s))
1306                         #
1307                         #   can we expect reduction of g?
1308                         #
1309                         # u_last=u_new
1310                         sub_iter+=1
1311                      self.relaxation=s
1312                #
1313                #  check for convergence:
1314                #
1315                norm_u_new = util.Lsup(u_new)
1316                p_new=p+self.relaxation*g_new
1317                norm_p_new = util.sqrt(self.inner(p_new,p_new))
1318                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))
1319    
1320                if self.iter>1:
1321                   dg2=g_new-g
1322                   df2=f_new-f
1323                   norm_dg2=util.sqrt(self.inner(dg2,dg2))
1324                   norm_df2=util.Lsup(df2)
1325                   # print norm_g_new, norm_g, norm_dg, norm_p, tolerance
1326                   tol_eq_g=tolerance*norm_dg2/(norm_g*abs(self.relaxation))*norm_p_new
1327                   tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new
1328                   if norm_g_new <= tol_eq_g and norm_f_new <= tol_eq_f:
1329                       converged=True
1330                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
1331            self.trace("convergence after %s steps."%self.iter)
1332            return u,p
1333    #   def solve(self,u0,p0,tolerance=1.e-6,iter_max=10,self.relaxation=1.):
1334    #      tol=1.e-2
1335    #      iter=0
1336    #      converged=False
1337    #      u=u0*1.
1338    #      p=p0*1.
1339    #      while not converged and iter<iter_max:
1340    #          du=self.solve_f(u,p,tol)
1341    #          u-=du
1342    #          norm_du=util.Lsup(du)
1343    #          norm_u=util.Lsup(u)
1344    #        
1345    #          dp=self.relaxation*self.solve_g(u,tol)
1346    #          p+=dp
1347    #          norm_dp=util.sqrt(self.inner(dp,dp))
1348    #          norm_p=util.sqrt(self.inner(p,p))
1349    #          print iter,"-th step rel. errror u,p= (%s,%s) (%s,%s)(%s,%s)"%(norm_du,norm_dp,norm_du/norm_u,norm_dp/norm_p,norm_u,norm_p)
1350    #          iter+=1
1351    #
1352    #          converged = (norm_du <= tolerance*norm_u) and  (norm_dp <= tolerance*norm_p)
1353    #      if converged:
1354    #          print "convergence after %s steps."%iter
1355    #      else:
1356    #          raise ArithmeticError("no convergence after %s steps."%iter)
1357    #
1358    #      return u,p
1359              
1360    def MaskFromBoundaryTag(function_space,*tags):
1361       """
1362       create a mask on the given function space which one for samples
1363       that touch the boundary tagged by tags.
1364    
1365       usage: m=MaskFromBoundaryTag(Solution(domain),"left", "right")
1366    
1367       @param function_space: a given function space
1368       @type function_space: L{escript.FunctionSpace}
1369       @param tags: boundray tags
1370       @type tags: C{str}
1371       @return: a mask which marks samples used by C{function_space} that are touching the
1372                boundary tagged by any of the given tags.
1373       @rtype: L{escript.Data} of rank 0
1374       """
1375       pde=linearPDEs.LinearPDE(function_space.getDomain(),numEquations=1, numSolutions=1)
1376       d=escript.Scalar(0.,escript.FunctionOnBoundary(function_space.getDomain()))
1377       for t in tags: d.setTaggedValue(t,1.)
1378       pde.setValue(y=d)
1379       out=util.whereNonZero(pde.getRightHandSide())
1380       if out.getFunctionSpace() == function_space:
1381          return out
1382       else:
1383          return util.whereNonZero(util.interpolate(out,function_space))
1384    
1385    
1386    

Legend:
Removed from v.782  
changed lines
  Added in v.1482

  ViewVC Help
Powered by ViewVC 1.1.26