/[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 1467 by gross, Wed Apr 2 08:10:37 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 PCG(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
481       """
482       Solver for
483    
484       M{Ax=b}
485    
486       with a symmetric and positive definite operator A (more details required!).
487       It uses the conjugate gradient method with preconditioner M providing an approximation of A.
488    
489       The iteration is terminated if the C{stoppingcriterium} function return C{True}.
490    
491       For details on the preconditioned conjugate gradient method see the book:
492    
493       Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
494       T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
495       C. Romine, and H. van der Vorst.
496    
497       @param b: the right hand side of the liner system. C{b} is altered.
498       @type b: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
499       @param Aprod: returns the value Ax
500       @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}.
501       @param Msolve: solves Mx=r
502       @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
503    type like argument C{x}.
504       @param bilinearform: inner product C{<x,r>}
505       @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}.
506       @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.
507       @type stoppingcriterium: function that returns C{True} or C{False}
508       @param x: an initial guess for the solution. If no C{x} is given 0*b is used.
509       @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
510       @param iter_max: maximum number of iteration steps.
511       @type iter_max: C{int}
512       @return: the solution approximation and the corresponding residual
513       @rtype: C{tuple}
514       @warning: C{b} and C{x} are altered.
515       """
516       iter=0
517       if x==None:
518          x=0*b
519       else:
520          b += (-1)*Aprod(x)
521       r=b
522       rhat=Msolve(r)
523       d = rhat
524       rhat_dot_r = bilinearform(rhat, r)
525       if rhat_dot_r<0: raise NegativeNorm,"negative norm."
526    
527       while not stoppingcriterium(math.sqrt(rhat_dot_r),r,x):
528           iter+=1
529           if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
530    
531           q=Aprod(d)
532           alpha = rhat_dot_r / bilinearform(d, q)
533           x += alpha * d
534           r += (-alpha) * q
535    
536           rhat=Msolve(r)
537           rhat_dot_r_new = bilinearform(rhat, r)
538           beta = rhat_dot_r_new / rhat_dot_r
539           rhat+=beta * d
540           d=rhat
541    
542           rhat_dot_r = rhat_dot_r_new
543           if rhat_dot_r<0: raise NegativeNorm,"negative norm."
544    
545       return x,r
546    
547    
548    ############################
549    # Added by Artak
550    #################################3
551    
552    #Apply a sequence of k Givens rotations, used within gmres codes
553    # vrot=givapp(c, s, vin, k)
554    def givapp(c,s,vin):
555        vrot=vin # warning: vin is altered!!!!
556        if isinstance(c,float):
557            vrot=[c*vrot[0]-s*vrot[1],s*vrot[0]+c*vrot[1]]
558        else:
559            for i in range(len(c)):
560                w1=c[i]*vrot[i]-s[i]*vrot[i+1]
561            w2=s[i]*vrot[i]+c[i]*vrot[i+1]
562                vrot[i:i+2]=w1,w2
563        return vrot
564    
565    def GMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
566       iter=0
567       r=Msolve(b)
568       r_dot_r = bilinearform(r, r)
569       if r_dot_r<0: raise NegativeNorm,"negative norm."
570       norm_b=math.sqrt(r_dot_r)
571    
572       if x==None:
573          x=0*b
574       else:
575          r=Msolve(b-Aprod(x))
576          r_dot_r = bilinearform(r, r)
577          if r_dot_r<0: raise NegativeNorm,"negative norm."
578      
579       h=numarray.zeros((iter_max,iter_max),numarray.Float64)
580       c=numarray.zeros(iter_max,numarray.Float64)
581       s=numarray.zeros(iter_max,numarray.Float64)
582       g=numarray.zeros(iter_max,numarray.Float64)
583       v=[]
584    
585       rho=math.sqrt(r_dot_r)
586       v.append(r/rho)
587       g[0]=rho
588    
589       while not stoppingcriterium(rho,norm_b):
590    
591        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
592    
593        
594        p=Msolve(Aprod(v[iter]))
595    
596        v.append(p)
597    
598        v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
599    
600    # Modified Gram-Schmidt
601        for j in range(iter+1):
602          h[j][iter]=bilinearform(v[j],v[iter+1])  
603          v[iter+1]+=(-1.)*h[j][iter]*v[j]
604          
605        h[iter+1][iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
606        v_norm2=h[iter+1][iter]
607    
608    
609    # Reorthogonalize if needed
610        if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
611         for j in range(iter+1):
612            hr=bilinearform(v[j],v[iter+1])
613                h[j][iter]=h[j][iter]+hr #vhat
614                v[iter+1] +=(-1.)*hr*v[j]
615    
616         v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
617         h[iter+1][iter]=v_norm2
618    
619    #   watch out for happy breakdown
620            if v_norm2 != 0:
621             v[iter+1]=v[iter+1]/h[iter+1][iter]
622    
623    #   Form and store the information for the new Givens rotation
624        if iter > 0 :
625            hhat=[]
626            for i in range(iter+1) : hhat.append(h[i][iter])
627            hhat=givapp(c[0:iter],s[0:iter],hhat);
628                for i in range(iter+1) : h[i][iter]=hhat[i]
629    
630        mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])
631        if mu!=0 :
632            c[iter]=h[iter][iter]/mu
633            s[iter]=-h[iter+1][iter]/mu
634            h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter]
635            h[iter+1][iter]=0.0
636            g[iter:iter+2]=givapp(c[iter],s[iter],g[iter:iter+2])
637    
638    # Update the residual norm
639            rho=abs(g[iter+1])
640        iter+=1
641    
642    # At this point either iter > iter_max or rho < tol.
643    # It's time to compute x and leave.        
644    
645       if iter > 0 :
646         y=numarray.zeros(iter,numarray.Float64)    
647         y[iter-1] = g[iter-1] / h[iter-1][iter-1]
648         if iter > 1 :  
649            i=iter-2  
650            while i>=0 :
651              y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]
652              i=i-1
653         xhat=v[iter-1]*y[iter-1]
654         for i in range(iter-1):
655        xhat += v[i]*y[i]
656       else : xhat=v[0]
657        
658       x += xhat
659    
660       return x
661        
662    #############################################
663    
664    class ArithmeticTuple(object):
665       """
666       tuple supporting inplace update x+=y and scaling x=a*y where x,y is an ArithmeticTuple and a is a float.
667    
668       example of usage:
669    
670       from esys.escript import Data
671       from numarray import array
672       a=Data(...)
673       b=array([1.,4.])
674       x=ArithmeticTuple(a,b)
675       y=5.*x
676    
677       """
678       def __init__(self,*args):
679           """
680           initialize object with elements args.
681    
682           @param args: tuple of object that support implace add (x+=y) and scaling (x=a*y)
683           """
684           self.__items=list(args)
685    
686       def __len__(self):
687           """
688           number of items
689    
690           @return: number of items
691           @rtype: C{int}
692           """
693           return len(self.__items)
694    
695       def __getitem__(self,index):
696           """
697           get an item
698    
699           @param index: item to be returned
700           @type index: C{int}
701           @return: item with index C{index}
702           """
703           return self.__items.__getitem__(index)
704    
705       def __mul__(self,other):
706           """
707           scaling from the right
708    
709           @param other: scaling factor
710           @type other: C{float}
711           @return: itemwise self*other
712           @rtype: L{ArithmeticTuple}
713           """
714           out=[]
715           for i in range(len(self)):
716               out.append(self[i]*other)
717           return ArithmeticTuple(*tuple(out))
718    
719       def __rmul__(self,other):
720           """
721           scaling from the left
722    
723           @param other: scaling factor
724           @type other: C{float}
725           @return: itemwise other*self
726           @rtype: L{ArithmeticTuple}
727           """
728           out=[]
729           for i in range(len(self)):
730               out.append(other*self[i])
731           return ArithmeticTuple(*tuple(out))
732    
733    #########################
734    # Added by Artak
735    #########################
736       def __div__(self,other):
737           """
738           dividing from the right
739    
740           @param other: scaling factor
741           @type other: C{float}
742           @return: itemwise self/other
743           @rtype: L{ArithmeticTuple}
744           """
745           out=[]
746           for i in range(len(self)):
747               out.append(self[i]/other)
748           return ArithmeticTuple(*tuple(out))
749    
750       def __rdiv__(self,other):
751           """
752           dividing from the left
753    
754           @param other: scaling factor
755           @type other: C{float}
756           @return: itemwise other/self
757           @rtype: L{ArithmeticTuple}
758           """
759           out=[]
760           for i in range(len(self)):
761               out.append(other/self[i])
762           return ArithmeticTuple(*tuple(out))
763      
764    ##########################################33
765    
766       def __iadd__(self,other):
767           """
768           in-place add of other to self
769    
770           @param other: increment
771           @type other: C{ArithmeticTuple}
772           """
773           if len(self) != len(other):
774               raise ValueError,"tuple length must match."
775           for i in range(len(self)):
776               self.__items[i]+=other[i]
777           return self
778    
779    class HomogeneousSaddlePointProblem(object):
780          """
781          This provides a framwork for solving homogeneous saddle point problem of the form
782    
783                 Av+B^*p=f
784                 Bv    =0
785    
786          for the unknowns v and p and given operators A and B and given right hand side f.
787          B^* is the adjoint operator of B is the given inner product.
788    
789          """
790          def __init__(self,**kwargs):
791            self.setTolerance()
792            self.setToleranceReductionFactor()
793    
794          def initialize(self):
795            """
796            initialize the problem (overwrite)
797            """
798            pass
799          def B(self,v):
800             """
801             returns Bv (overwrite)
802             @rtype: equal to the type of p
803    
804             @note: boundary conditions on p should be zero!
805             """
806             pass
807    
808          def inner(self,p0,p1):
809             """
810             returns inner product of two element p0 and p1  (overwrite)
811            
812             @type p0: equal to the type of p
813             @type p1: equal to the type of p
814             @rtype: C{float}
815    
816             @rtype: equal to the type of p
817             """
818             pass
819    
820          def solve_A(self,u,p):
821             """
822             solves Av=f-Au-B^*p with accuracy self.getReducedTolerance() (overwrite)
823    
824             @rtype: equal to the type of v
825             @note: boundary conditions on v should be zero!
826             """
827             pass
828    
829          def solve_prec(self,p):
830             """
831             provides a preconditioner for BA^{-1}B^* with accuracy self.getReducedTolerance() (overwrite)
832    
833             @rtype: equal to the type of p
834             """
835             pass
836    
837          def stoppingcriterium(self,Bv,v,p):
838             """
839             returns a True if iteration is terminated. (overwrite)
840    
841             @rtype: C{bool}
842             """
843             pass
844                
845          def __inner(self,p,r):
846             return self.inner(p,r[1])
847    
848          def __inner_p(self,p1,p2):
849             return self.inner(p1,p2)
850    
851          def __stoppingcriterium(self,norm_r,r,p):
852              return self.stoppingcriterium(r[1],r[0],p)
853    
854          def __stoppingcriterium_GMRES(self,norm_r,norm_b):
855              return self.stoppingcriterium_GMRES(norm_r,norm_b)
856    
857          def setTolerance(self,tolerance=1.e-8):
858                  self.__tol=tolerance
859          def getTolerance(self):
860                  return self.__tol
861          def setToleranceReductionFactor(self,reduction=0.01):
862                  self.__reduction=reduction
863          def getSubProblemTolerance(self):
864                  return self.__reduction*self.getTolerance()
865    
866          def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='GMRES'):
867                  """
868                  solves the saddle point problem using initial guesses v and p.
869    
870                  @param max_iter: maximum number of iteration steps.
871                  """
872                  self.verbose=verbose
873                  self.show_details=show_details and self.verbose
874    
875              # Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask)      
876    
877              
878              self.__z=v+self.solve_A(v,p*0)
879    
880                  Bz=self.B(self.__z)
881                  #
882              #   solve BA^-1B^*p = Bz
883                  #
884                  #   note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv
885                  #
886                  #   with                    Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask)
887                  #                           A(v-z)=Az-B^*p-Az = f -Az - B^*p (v-z=0 on fixed_u_mask)
888                  #
889                  self.iter=0
890              if solver=='GMRES':      
891                    if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter
892                    p=GMRES(Bz,self.__Aprod_GMRES,self.__Msolve_GMRES,self.__inner_p,self.__stoppingcriterium_GMRES,iter_max=max_iter, x=p*1.)
893                    # solve Au=f-B^*p
894                    #       A(u-v)=f-B^*p-Av
895                    #       u=v+(u-v)
896            u=v+self.solve_A(v,p)
897        
898                  else:
899                    if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter
900                    p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p)
901                u=r[0]  
902    
903              return u,p
904    
905          def __Msolve(self,r):
906              return self.solve_prec(r[1])
907    
908          def __Msolve_GMRES(self,r):
909              return self.solve_prec(r)
910    
911    
912          def __Aprod(self,p):
913              # return BA^-1B*p
914              #solve Av =-B^*p as Av =f-Az-B^*p
915              v=self.solve_A(self.__z,p)
916              return ArithmeticTuple(v, -self.B(v))
917    
918          def __Aprod_GMRES(self,p):
919              # return BA^-1B*p
920              #solve Av =-B^*p as Av =f-Az-B^*p
921          v=self.solve_A(self.__z,p)
922              return -self.B(v)
923    
924  class SaddlePointProblem(object):  class SaddlePointProblem(object):
925     """     """
926     This implements a solver for a saddlepoint problem     This implements a solver for a saddlepoint problem
# Line 363  class SaddlePointProblem(object): Line 930  class SaddlePointProblem(object):
930    
931     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:
932    
933     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})}
934     M{Q_g (p^{k+1}-p^{k}) =   g(u^{k+1})}     M{Q_g (p^{k+1}-p^{k}) =   g(u^{k+1})}
935    
936     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 942  class SaddlePointProblem(object):
942         """         """
943         initializes the problem         initializes the problem
944    
945         @parm verbose: switches on the printing out some information         @param verbose: switches on the printing out some information
946         @type verbose: C{bool}         @type verbose: C{bool}
947         @note: this method may be overwritten by a particular saddle point problem         @note: this method may be overwritten by a particular saddle point problem
948         """         """
949           if not isinstance(verbose,bool):
950                raise TypeError("verbose needs to be of type bool.")
951         self.__verbose=verbose         self.__verbose=verbose
952         self.relaxation=1.         self.relaxation=1.
953    
# Line 386  class SaddlePointProblem(object): Line 955  class SaddlePointProblem(object):
955         """         """
956         prints text if verbose has been set         prints text if verbose has been set
957    
958         @parm text: a text message         @param text: a text message
959         @type text: C{str}         @type text: C{str}
960         """         """
961         if self.__verbose: print "%s: %s"%(str(self),text)         if self.__verbose: print "%s: %s"%(str(self),text)
# Line 539  class SaddlePointProblem(object): Line 1108  class SaddlePointProblem(object):
1108              norm_u_new = util.Lsup(u_new)              norm_u_new = util.Lsup(u_new)
1109              p_new=p+self.relaxation*g_new              p_new=p+self.relaxation*g_new
1110              norm_p_new = util.sqrt(self.inner(p_new,p_new))              norm_p_new = util.sqrt(self.inner(p_new,p_new))
1111              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))
1112    
1113              if self.iter>1:              if self.iter>1:
1114                 dg2=g_new-g                 dg2=g_new-g
# Line 551  class SaddlePointProblem(object): Line 1120  class SaddlePointProblem(object):
1120                 tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new                 tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new
1121                 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:
1122                     converged=True                     converged=True
                    break  
1123              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
1124          self.trace("convergence after %s steps."%self.iter)          self.trace("convergence after %s steps."%self.iter)
1125          return u,p          return u,p
# Line 582  class SaddlePointProblem(object): Line 1150  class SaddlePointProblem(object):
1150  #  #
1151  #      return u,p  #      return u,p
1152                        
1153  # vim: expandtab shiftwidth=4:  def MaskFromBoundaryTag(function_space,*tags):
1154       """
1155       create a mask on the given function space which one for samples
1156       that touch the boundary tagged by tags.
1157    
1158       usage: m=MaskFromBoundaryTag(Solution(domain),"left", "right")
1159    
1160       @param function_space: a given function space
1161       @type function_space: L{escript.FunctionSpace}
1162       @param tags: boundray tags
1163       @type tags: C{str}
1164       @return: a mask which marks samples used by C{function_space} that are touching the
1165                boundary tagged by any of the given tags.
1166       @rtype: L{escript.Data} of rank 0
1167       """
1168       pde=linearPDEs.LinearPDE(function_space.getDomain(),numEquations=1, numSolutions=1)
1169       d=escript.Scalar(0.,escript.FunctionOnBoundary(function_space.getDomain()))
1170       for t in tags: d.setTaggedValue(t,1.)
1171       pde.setValue(y=d)
1172       out=util.whereNonZero(pde.getRightHandSide())
1173       if out.getFunctionSpace() == function_space:
1174          return out
1175       else:
1176          return util.whereNonZero(util.interpolate(out,function_space))
1177    
1178    
1179    

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

  ViewVC Help
Powered by ViewVC 1.1.26