/[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 1469 by gross, Thu Apr 3 05:16:56 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):
583       iter=0
584       r=Msolve(b)
585       r_dot_r = bilinearform(r, r)
586       if r_dot_r<0: raise NegativeNorm,"negative norm."
587       norm_b=math.sqrt(r_dot_r)
588    
589       if x==None:
590          x=0*b
591       else:
592          r=Msolve(b-Aprod(x))
593          r_dot_r = bilinearform(r, r)
594          if r_dot_r<0: raise NegativeNorm,"negative norm."
595      
596       h=numarray.zeros((iter_max,iter_max),numarray.Float64)
597       c=numarray.zeros(iter_max,numarray.Float64)
598       s=numarray.zeros(iter_max,numarray.Float64)
599       g=numarray.zeros(iter_max,numarray.Float64)
600       v=[]
601    
602       rho=math.sqrt(r_dot_r)
603       v.append(r/rho)
604       g[0]=rho
605    
606       while not stoppingcriterium(rho,norm_b):
607    
608        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
609    
610        
611        p=Msolve(Aprod(v[iter]))
612    
613        v.append(p)
614    
615        v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
616    
617    # Modified Gram-Schmidt
618        for j in range(iter+1):
619          h[j][iter]=bilinearform(v[j],v[iter+1])  
620          v[iter+1]+=(-1.)*h[j][iter]*v[j]
621          
622        h[iter+1][iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
623        v_norm2=h[iter+1][iter]
624    
625    
626    # Reorthogonalize if needed
627        if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
628         for j in range(iter+1):
629            hr=bilinearform(v[j],v[iter+1])
630                h[j][iter]=h[j][iter]+hr #vhat
631                v[iter+1] +=(-1.)*hr*v[j]
632    
633         v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
634         h[iter+1][iter]=v_norm2
635    
636    #   watch out for happy breakdown
637            if v_norm2 != 0:
638             v[iter+1]=v[iter+1]/h[iter+1][iter]
639    
640    #   Form and store the information for the new Givens rotation
641        if iter > 0 :
642            hhat=[]
643            for i in range(iter+1) : hhat.append(h[i][iter])
644            hhat=givapp(c[0:iter],s[0:iter],hhat);
645                for i in range(iter+1) : h[i][iter]=hhat[i]
646    
647        mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])
648        if mu!=0 :
649            c[iter]=h[iter][iter]/mu
650            s[iter]=-h[iter+1][iter]/mu
651            h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter]
652            h[iter+1][iter]=0.0
653            g[iter:iter+2]=givapp(c[iter],s[iter],g[iter:iter+2])
654    
655    # Update the residual norm
656            rho=abs(g[iter+1])
657        iter+=1
658    
659    # At this point either iter > iter_max or rho < tol.
660    # It's time to compute x and leave.        
661    
662       if iter > 0 :
663         y=numarray.zeros(iter,numarray.Float64)    
664         y[iter-1] = g[iter-1] / h[iter-1][iter-1]
665         if iter > 1 :  
666            i=iter-2  
667            while i>=0 :
668              y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]
669              i=i-1
670         xhat=v[iter-1]*y[iter-1]
671         for i in range(iter-1):
672        xhat += v[i]*y[i]
673       else : xhat=v[0]
674        
675       x += xhat
676    
677       return x
678        
679    #############################################
680    
681    class ArithmeticTuple(object):
682       """
683       tuple supporting inplace update x+=y and scaling x=a*y where x,y is an ArithmeticTuple and a is a float.
684    
685       example of usage:
686    
687       from esys.escript import Data
688       from numarray import array
689       a=Data(...)
690       b=array([1.,4.])
691       x=ArithmeticTuple(a,b)
692       y=5.*x
693    
694       """
695       def __init__(self,*args):
696           """
697           initialize object with elements args.
698    
699           @param args: tuple of object that support implace add (x+=y) and scaling (x=a*y)
700           """
701           self.__items=list(args)
702    
703       def __len__(self):
704           """
705           number of items
706    
707           @return: number of items
708           @rtype: C{int}
709           """
710           return len(self.__items)
711    
712       def __getitem__(self,index):
713           """
714           get an item
715    
716           @param index: item to be returned
717           @type index: C{int}
718           @return: item with index C{index}
719           """
720           return self.__items.__getitem__(index)
721    
722       def __mul__(self,other):
723           """
724           scaling from the right
725    
726           @param other: scaling factor
727           @type other: C{float}
728           @return: itemwise self*other
729           @rtype: L{ArithmeticTuple}
730           """
731           out=[]
732           for i in range(len(self)):
733               out.append(self[i]*other)
734           return ArithmeticTuple(*tuple(out))
735    
736       def __rmul__(self,other):
737           """
738           scaling from the left
739    
740           @param other: scaling factor
741           @type other: C{float}
742           @return: itemwise other*self
743           @rtype: L{ArithmeticTuple}
744           """
745           out=[]
746           for i in range(len(self)):
747               out.append(other*self[i])
748           return ArithmeticTuple(*tuple(out))
749    
750    #########################
751    # Added by Artak
752    #########################
753       def __div__(self,other):
754           """
755           dividing from the right
756    
757           @param other: scaling factor
758           @type other: C{float}
759           @return: itemwise self/other
760           @rtype: L{ArithmeticTuple}
761           """
762           out=[]
763           for i in range(len(self)):
764               out.append(self[i]/other)
765           return ArithmeticTuple(*tuple(out))
766    
767       def __rdiv__(self,other):
768           """
769           dividing from the left
770    
771           @param other: scaling factor
772           @type other: C{float}
773           @return: itemwise other/self
774           @rtype: L{ArithmeticTuple}
775           """
776           out=[]
777           for i in range(len(self)):
778               out.append(other/self[i])
779           return ArithmeticTuple(*tuple(out))
780      
781    ##########################################33
782    
783       def __iadd__(self,other):
784           """
785           in-place add of other to self
786    
787           @param other: increment
788           @type other: C{ArithmeticTuple}
789           """
790           if len(self) != len(other):
791               raise ValueError,"tuple length must match."
792           for i in range(len(self)):
793               self.__items[i]+=other[i]
794           return self
795    
796    class HomogeneousSaddlePointProblem(object):
797          """
798          This provides a framwork for solving homogeneous saddle point problem of the form
799    
800                 Av+B^*p=f
801                 Bv    =0
802    
803          for the unknowns v and p and given operators A and B and given right hand side f.
804          B^* is the adjoint operator of B is the given inner product.
805    
806          """
807          def __init__(self,**kwargs):
808            self.setTolerance()
809            self.setToleranceReductionFactor()
810    
811          def initialize(self):
812            """
813            initialize the problem (overwrite)
814            """
815            pass
816          def B(self,v):
817             """
818             returns Bv (overwrite)
819             @rtype: equal to the type of p
820    
821             @note: boundary conditions on p should be zero!
822             """
823             pass
824    
825          def inner(self,p0,p1):
826             """
827             returns inner product of two element p0 and p1  (overwrite)
828            
829             @type p0: equal to the type of p
830             @type p1: equal to the type of p
831             @rtype: C{float}
832    
833             @rtype: equal to the type of p
834             """
835             pass
836    
837          def solve_A(self,u,p):
838             """
839             solves Av=f-Au-B^*p with accuracy self.getReducedTolerance() (overwrite)
840    
841             @rtype: equal to the type of v
842             @note: boundary conditions on v should be zero!
843             """
844             pass
845    
846          def solve_prec(self,p):
847             """
848             provides a preconditioner for BA^{-1}B^* with accuracy self.getReducedTolerance() (overwrite)
849    
850             @rtype: equal to the type of p
851             """
852             pass
853    
854          def stoppingcriterium(self,Bv,v,p):
855             """
856             returns a True if iteration is terminated. (overwrite)
857    
858             @rtype: C{bool}
859             """
860             pass
861                
862          def __inner(self,p,r):
863             return self.inner(p,r[1])
864    
865          def __inner_p(self,p1,p2):
866             return self.inner(p1,p2)
867    
868          def __stoppingcriterium(self,norm_r,r,p):
869              return self.stoppingcriterium(r[1],r[0],p)
870    
871          def __stoppingcriterium_GMRES(self,norm_r,norm_b):
872              return self.stoppingcriterium_GMRES(norm_r,norm_b)
873    
874          def setTolerance(self,tolerance=1.e-8):
875                  self.__tol=tolerance
876          def getTolerance(self):
877                  return self.__tol
878          def setToleranceReductionFactor(self,reduction=0.01):
879                  self.__reduction=reduction
880          def getSubProblemTolerance(self):
881                  return self.__reduction*self.getTolerance()
882    
883          def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='GMRES'):
884                  """
885                  solves the saddle point problem using initial guesses v and p.
886    
887                  @param max_iter: maximum number of iteration steps.
888                  """
889                  self.verbose=verbose
890                  self.show_details=show_details and self.verbose
891    
892                  # assume p is known: then v=A^-1(f-B^*p)
893                  # which leads to BA^-1B^*p = BA^-1f  
894    
895              # Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask)      
896    
897              
898              self.__z=v+self.solve_A(v,p*0)
899    
900                  Bz=self.B(self.__z)
901                  #
902              #   solve BA^-1B^*p = Bz
903                  #
904                  #   note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv
905                  #
906                  #   with                    Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask)
907                  #                           A(v-z)=Az-B^*p-Az = f -Az - B^*p (v-z=0 on fixed_u_mask)
908                  #
909                  self.iter=0
910              if solver=='GMRES':      
911                    if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter
912                    p=GMRES(Bz,self.__Aprod_GMRES,self.__Msolve_GMRES,self.__inner_p,self.__stoppingcriterium_GMRES,iter_max=max_iter, x=p*1.)
913                    # solve Au=f-B^*p
914                    #       A(u-v)=f-B^*p-Av
915                    #       u=v+(u-v)
916            u=v+self.solve_A(v,p)
917        
918                  else:
919                    if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter
920                    p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p)
921                u=r[0]  
922                  print "div(u)=",util.Lsup(self.B(u)),util.Lsup(u)
923    
924              return u,p
925    
926          def __Msolve(self,r):
927              return self.solve_prec(r[1])
928    
929          def __Msolve_GMRES(self,r):
930              return self.solve_prec(r)
931    
932    
933          def __Aprod(self,p):
934              # return BA^-1B*p
935              #solve Av =-B^*p as Av =f-Az-B^*p
936              v=self.solve_A(self.__z,-p)
937              return ArithmeticTuple(v, self.B(v))
938    
939          def __Aprod_GMRES(self,p):
940              # return BA^-1B*p
941              #solve Av =-B^*p as Av =f-Az-B^*p
942          v=self.solve_A(self.__z,-p)
943              return self.B(v)
944    
945  class SaddlePointProblem(object):  class SaddlePointProblem(object):
946     """     """
947     This implements a solver for a saddlepoint problem     This implements a solver for a saddlepoint problem
# Line 363  class SaddlePointProblem(object): Line 951  class SaddlePointProblem(object):
951    
952     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:
953    
954     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})}
955     M{Q_g (p^{k+1}-p^{k}) =   g(u^{k+1})}     M{Q_g (p^{k+1}-p^{k}) =   g(u^{k+1})}
956    
957     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 963  class SaddlePointProblem(object):
963         """         """
964         initializes the problem         initializes the problem
965    
966         @parm verbose: switches on the printing out some information         @param verbose: switches on the printing out some information
967         @type verbose: C{bool}         @type verbose: C{bool}
968         @note: this method may be overwritten by a particular saddle point problem         @note: this method may be overwritten by a particular saddle point problem
969         """         """
970           if not isinstance(verbose,bool):
971                raise TypeError("verbose needs to be of type bool.")
972         self.__verbose=verbose         self.__verbose=verbose
973         self.relaxation=1.         self.relaxation=1.
974    
# Line 386  class SaddlePointProblem(object): Line 976  class SaddlePointProblem(object):
976         """         """
977         prints text if verbose has been set         prints text if verbose has been set
978    
979         @parm text: a text message         @param text: a text message
980         @type text: C{str}         @type text: C{str}
981         """         """
982         if self.__verbose: print "%s: %s"%(str(self),text)         if self.__verbose: print "%s: %s"%(str(self),text)
# Line 539  class SaddlePointProblem(object): Line 1129  class SaddlePointProblem(object):
1129              norm_u_new = util.Lsup(u_new)              norm_u_new = util.Lsup(u_new)
1130              p_new=p+self.relaxation*g_new              p_new=p+self.relaxation*g_new
1131              norm_p_new = util.sqrt(self.inner(p_new,p_new))              norm_p_new = util.sqrt(self.inner(p_new,p_new))
1132              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))
1133    
1134              if self.iter>1:              if self.iter>1:
1135                 dg2=g_new-g                 dg2=g_new-g
# Line 551  class SaddlePointProblem(object): Line 1141  class SaddlePointProblem(object):
1141                 tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new                 tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new
1142                 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:
1143                     converged=True                     converged=True
                    break  
1144              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
1145          self.trace("convergence after %s steps."%self.iter)          self.trace("convergence after %s steps."%self.iter)
1146          return u,p          return u,p
# Line 582  class SaddlePointProblem(object): Line 1171  class SaddlePointProblem(object):
1171  #  #
1172  #      return u,p  #      return u,p
1173                        
1174  # vim: expandtab shiftwidth=4:  def MaskFromBoundaryTag(function_space,*tags):
1175       """
1176       create a mask on the given function space which one for samples
1177       that touch the boundary tagged by tags.
1178    
1179       usage: m=MaskFromBoundaryTag(Solution(domain),"left", "right")
1180    
1181       @param function_space: a given function space
1182       @type function_space: L{escript.FunctionSpace}
1183       @param tags: boundray tags
1184       @type tags: C{str}
1185       @return: a mask which marks samples used by C{function_space} that are touching the
1186                boundary tagged by any of the given tags.
1187       @rtype: L{escript.Data} of rank 0
1188       """
1189       pde=linearPDEs.LinearPDE(function_space.getDomain(),numEquations=1, numSolutions=1)
1190       d=escript.Scalar(0.,escript.FunctionOnBoundary(function_space.getDomain()))
1191       for t in tags: d.setTaggedValue(t,1.)
1192       pde.setValue(y=d)
1193       out=util.whereNonZero(pde.getRightHandSide())
1194       if out.getFunctionSpace() == function_space:
1195          return out
1196       else:
1197          return util.whereNonZero(util.interpolate(out,function_space))
1198    
1199    
1200    

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

  ViewVC Help
Powered by ViewVC 1.1.26