/[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 1414 by gross, Thu Feb 14 10:01:43 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  class TimeIntegrationManager:  class TimeIntegrationManager:
52    """    """
# Line 131  class Projector: Line 147  class Projector:
147      self.__pde.setValue(D = 1.)      self.__pde.setValue(D = 1.)
148      return      return
149    
   def __del__(self):  
     return  
   
150    def __call__(self, input_data):    def __call__(self, input_data):
151      """      """
152      Projects input_data onto a continuous function      Projects input_data onto a continuous function
# Line 141  class Projector: Line 154  class Projector:
154      @param input_data: The input_data to be projected.      @param input_data: The input_data to be projected.
155      """      """
156      out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())      out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())
157        self.__pde.setValue(Y = escript.Data(), Y_reduced = escript.Data())
158      if input_data.getRank()==0:      if input_data.getRank()==0:
159          self.__pde.setValue(Y = input_data)          self.__pde.setValue(Y = input_data)
160          out=self.__pde.getSolution()          out=self.__pde.getSolution()
# Line 296  class Locator: Line 310  class Locator:
310         Initializes a Locator to access values in Data objects on the Doamin         Initializes a Locator to access values in Data objects on the Doamin
311         or FunctionSpace where for the sample point which         or FunctionSpace where for the sample point which
312         closest to the given point x.         closest to the given point x.
313    
314           @param where: function space
315           @type where: L{escript.FunctionSpace}
316           @param x: coefficient of the solution.
317           @type x: L{numarray.NumArray} or C{list} of L{numarray.NumArray}
318         """         """
319         if isinstance(where,escript.FunctionSpace):         if isinstance(where,escript.FunctionSpace):
320            self.__function_space=where            self.__function_space=where
321         else:         else:
322            self.__function_space=escript.ContinuousFunction(where)            self.__function_space=escript.ContinuousFunction(where)
323         self.__id=util.length(self.__function_space.getX()-x[:self.__function_space.getDim()]).mindp()         if isinstance(x, list):
324               self.__id=[]
325               for p in x:
326                  self.__id.append(util.length(self.__function_space.getX()-p[:self.__function_space.getDim()]).minGlobalDataPoint())
327           else:
328               self.__id=util.length(self.__function_space.getX()-x[:self.__function_space.getDim()]).minGlobalDataPoint()
329    
330       def __str__(self):       def __str__(self):
331         """         """
332         Returns the coordinates of the Locator as a string.         Returns the coordinates of the Locator as a string.
333         """         """
334         return "<Locator %s>"%str(self.getX())         x=self.getX()
335           if instance(x,list):
336              out="["
337              first=True
338              for xx in x:
339                if not first:
340                    out+=","
341                else:
342                    first=False
343                out+=str(xx)
344              out+="]>"
345           else:
346              out=str(x)
347           return out
348    
349         def getX(self):
350            """
351        Returns the exact coordinates of the Locator.
352        """
353            return self(self.getFunctionSpace().getX())
354    
355       def getFunctionSpace(self):       def getFunctionSpace(self):
356          """          """
# Line 315  class Locator: Line 358  class Locator:
358      """      """
359          return self.__function_space          return self.__function_space
360    
361       def getId(self):       def getId(self,item=None):
362          """          """
363      Returns the identifier of the location.      Returns the identifier of the location.
364      """      """
365          return self.__id          if item == None:
366               return self.__id
367            else:
368               if isinstance(self.__id,list):
369                  return self.__id[item]
370               else:
371                  return self.__id
372    
      def getX(self):  
         """  
     Returns the exact coordinates of the Locator.  
     """  
         return self(self.getFunctionSpace().getX())  
373    
374       def __call__(self,data):       def __call__(self,data):
375          """          """
# Line 341  class Locator: Line 385  class Locator:
385      """      """
386          if isinstance(data,escript.Data):          if isinstance(data,escript.Data):
387             if data.getFunctionSpace()==self.getFunctionSpace():             if data.getFunctionSpace()==self.getFunctionSpace():
388               out=data.convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])               dat=data
              #out=data.convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1],self.getId()[2])  
389             else:             else:
390               out=data.interpolate(self.getFunctionSpace()).convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])               dat=data.interpolate(self.getFunctionSpace())
391               #out=data.interpolate(self.getFunctionSpace()).convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1],self.getId()[2])             id=self.getId()
392             if data.getRank()==0:             r=data.getRank()
393                return out[0]             if isinstance(id,list):
394                   out=[]
395                   for i in id:
396                      o=data.getValueOfGlobalDataPoint(*i)
397                      if data.getRank()==0:
398                         out.append(o[0])
399                      else:
400                         out.append(o)
401                   return out
402             else:             else:
403                return out               out=data.getValueOfGlobalDataPoint(*id)
404                 if data.getRank()==0:
405                    return out[0]
406                 else:
407                    return out
408          else:          else:
409             return data             return data
410    
411  # vim: expandtab shiftwidth=4:  class SolverSchemeException(Exception):
412       """
413       exceptions thrown by solvers
414       """
415       pass
416    
417    class IndefinitePreconditioner(SolverSchemeException):
418       """
419       the preconditioner is not positive definite.
420       """
421       pass
422    class MaxIterReached(SolverSchemeException):
423       """
424       maxium number of iteration steps is reached.
425       """
426       pass
427    class IterationBreakDown(SolverSchemeException):
428       """
429       iteration scheme econouters an incurable breakdown.
430       """
431       pass
432    class NegativeNorm(SolverSchemeException):
433       """
434       a norm calculation returns a negative norm.
435       """
436       pass
437    
438    class IterationHistory(object):
439       """
440       The IterationHistory class is used to define a stopping criterium. It keeps track of the
441       residual norms. The stoppingcriterium indicates termination if the residual norm has been reduced by
442       a given tolerance.
443       """
444       def __init__(self,tolerance=math.sqrt(util.EPSILON),verbose=False):
445          """
446          Initialization
447    
448          @param tolerance: tolerance
449          @type tolerance: positive C{float}
450          @param verbose: switches on the printing out some information
451          @type verbose: C{bool}
452          """
453          if not tolerance>0.:
454              raise ValueError,"tolerance needs to be positive."
455          self.tolerance=tolerance
456          self.verbose=verbose
457          self.history=[]
458       def stoppingcriterium(self,norm_r,r,x):
459           """
460           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.
461    
462          
463           @param norm_r: current residual norm
464           @type norm_r: non-negative C{float}
465           @param r: current residual (not used)
466           @param x: current solution approximation (not used)
467           @return: C{True} is the stopping criterium is fullfilled. Otherwise C{False} is returned.
468           @rtype: C{bool}
469    
470           """
471           self.history.append(norm_r)
472           if self.verbose: print "iter: %s:  inner(rhat,r) = %e"%(len(self.history)-1, self.history[-1])
473           return self.history[-1]<=self.tolerance * self.history[0]
474    
475    def PCG(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
476       """
477       Solver for
478    
479       M{Ax=b}
480    
481       with a symmetric and positive definite operator A (more details required!).
482       It uses the conjugate gradient method with preconditioner M providing an approximation of A.
483    
484       The iteration is terminated if the C{stoppingcriterium} function return C{True}.
485    
486       For details on the preconditioned conjugate gradient method see the book:
487    
488       Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
489       T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
490       C. Romine, and H. van der Vorst.
491    
492       @param b: the right hand side of the liner system. C{b} is altered.
493       @type b: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
494       @param Aprod: returns the value Ax
495       @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}.
496       @param Msolve: solves Mx=r
497       @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
498    type like argument C{x}.
499       @param bilinearform: inner product C{<x,r>}
500       @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}.
501       @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.
502       @type stoppingcriterium: function that returns C{True} or C{False}
503       @param x: an initial guess for the solution. If no C{x} is given 0*b is used.
504       @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
505       @param iter_max: maximum number of iteration steps.
506       @type iter_max: C{int}
507       @return: the solution approximation and the corresponding residual
508       @rtype: C{tuple}
509       @warning: C{b} and C{x} are altered.
510       """
511       iter=0
512       if x==None:
513          x=0*b
514       else:
515          b += (-1)*Aprod(x)
516       r=b
517       rhat=Msolve(r)
518       d = rhat
519       rhat_dot_r = bilinearform(rhat, r)
520       if rhat_dot_r<0: raise NegativeNorm,"negative norm."
521    
522       while not stoppingcriterium(math.sqrt(rhat_dot_r),r,x):
523           iter+=1
524           if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
525    
526           q=Aprod(d)
527           alpha = rhat_dot_r / bilinearform(d, q)
528           x += alpha * d
529           r += (-alpha) * q
530    
531           rhat=Msolve(r)
532           rhat_dot_r_new = bilinearform(rhat, r)
533           beta = rhat_dot_r_new / rhat_dot_r
534           rhat+=beta * d
535           d=rhat
536    
537           rhat_dot_r = rhat_dot_r_new
538           if rhat_dot_r<0: raise NegativeNorm,"negative norm."
539    
540       return x,r
541    
542    class ArithmeticTuple(object):
543       """
544       tuple supporting inplace update x+=y and scaling x=a*y where x,y is an ArithmeticTuple and a is a float.
545    
546       example of usage:
547    
548       from esys.escript import Data
549       from numarray import array
550       a=Data(...)
551       b=array([1.,4.])
552       x=ArithmeticTuple(a,b)
553       y=5.*x
554    
555       """
556       def __init__(self,*args):
557           """
558           initialize object with elements args.
559    
560           @param args: tuple of object that support implace add (x+=y) and scaling (x=a*y)
561           """
562           self.__items=list(args)
563    
564       def __len__(self):
565           """
566           number of items
567    
568           @return: number of items
569           @rtype: C{int}
570           """
571           return len(self.__items)
572    
573       def __getitem__(self,index):
574           """
575           get an item
576    
577           @param index: item to be returned
578           @type index: C{int}
579           @return: item with index C{index}
580           """
581           return self.__items.__getitem__(index)
582    
583       def __mul__(self,other):
584           """
585           scaling from the right
586    
587           @param other: scaling factor
588           @type other: C{float}
589           @return: itemwise self*other
590           @rtype: L{ArithmeticTuple}
591           """
592           out=[]
593           for i in range(len(self)):
594               out.append(self[i]*other)
595           return ArithmeticTuple(*tuple(out))
596    
597       def __rmul__(self,other):
598           """
599           scaling from the left
600    
601           @param other: scaling factor
602           @type other: C{float}
603           @return: itemwise other*self
604           @rtype: L{ArithmeticTuple}
605           """
606           out=[]
607           for i in range(len(self)):
608               out.append(other*self[i])
609           return ArithmeticTuple(*tuple(out))
610    
611       def __iadd__(self,other):
612           """
613           in-place add of other to self
614    
615           @param other: increment
616           @type other: C{ArithmeticTuple}
617           """
618           if len(self) != len(other):
619               raise ValueError,"tuple length must match."
620           for i in range(len(self)):
621               self.__items[i]+=other[i]
622           return self
623    
624    class HomogeneousSaddlePointProblem(object):
625          """
626          This provides a framwork for solving homogeneous saddle point problem of the form
627    
628                 Av+B^*p=f
629                 Bv    =0
630    
631          for the unknowns v and p and given operators A and B and given right hand side f.
632          B^* is the adjoint operator of B is the given inner product.
633    
634          """
635          def __init__(self,**kwargs):
636            self.setTolerance()
637            self.setToleranceReductionFactor()
638    
639          def initialize(self):
640            """
641            initialize the problem (overwrite)
642            """
643            pass
644          def B(self,v):
645             """
646             returns Bv (overwrite)
647             @rtype: equal to the type of p
648    
649             @note: boundary conditions on p should be zero!
650             """
651             pass
652    
653          def inner(self,p0,p1):
654             """
655             returns inner product of two element p0 and p1  (overwrite)
656            
657             @type p0: equal to the type of p
658             @type p1: equal to the type of p
659             @rtype: C{float}
660    
661             @rtype: equal to the type of p
662             """
663             pass
664    
665          def solve_A(self,u,p):
666             """
667             solves Av=f-Au-B^*p with accuracy self.getReducedTolerance() (overwrite)
668    
669             @rtype: equal to the type of v
670             @note: boundary conditions on v should be zero!
671             """
672             pass
673    
674          def solve_prec(self,p):
675             """
676             provides a preconditioner for BA^{-1}B^* with accuracy self.getReducedTolerance() (overwrite)
677    
678             @rtype: equal to the type of p
679             """
680             pass
681    
682          def stoppingcriterium(self,Bv,v,p):
683             """
684             returns a True if iteration is terminated. (overwrite)
685    
686             @rtype: C{bool}
687             """
688             pass
689                
690          def __inner(self,p,r):
691             return self.inner(p,r[1])
692    
693          def __stoppingcriterium(self,norm_r,r,p):
694              return self.stoppingcriterium(r[1],r[0],p)
695    
696          def setTolerance(self,tolerance=1.e-8):
697                  self.__tol=tolerance
698          def getTolerance(self):
699                  return self.__tol
700          def setToleranceReductionFactor(self,reduction=0.01):
701                  self.__reduction=reduction
702          def getSubProblemTolerance(self):
703                  return self.__reduction*self.getTolerance()
704    
705          def solve(self,v,p,max_iter=20, verbose=False, show_details=False):
706                  """
707                  solves the saddle point problem using initial guesses v and p.
708    
709                  @param max_iter: maximum number of iteration steps.
710                  """
711                  self.verbose=verbose
712                  self.show_details=show_details and self.verbose
713    
714              # Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask)      
715    
716              self.__z=v+self.solve_A(v,p*0)
717                  Bz=self.B(self.__z)
718                  #
719              #   solve BA^-1B^*p = Bz
720                  #
721                  #   note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv
722                  #
723                  #   with                    Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask)
724                  #                           A(v-z)=Az-B^*p-Az = f -Az - B^*p (v-z=0 on fixed_u_mask)
725                  #
726                  self.iter=0
727                  if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter
728                  p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p)
729              return r[0],p
730    
731          def __Msolve(self,r):
732              return self.solve_prec(r[1])
733    
734          def __Aprod(self,p):
735              # return BA^-1B*p
736              #solve Av =-B^*p as Av =f-Az-B^*p
737              v=self.solve_A(self.__z,p)
738              return ArithmeticTuple(v, self.B(v))
739    
740    
741    class SaddlePointProblem(object):
742       """
743       This implements a solver for a saddlepoint problem
744    
745       M{f(u,p)=0}
746       M{g(u)=0}
747    
748       for u and p. The problem is solved with an inexact Uszawa scheme for p:
749    
750       M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})}
751       M{Q_g (p^{k+1}-p^{k}) =   g(u^{k+1})}
752    
753       where Q_f is an approximation of the Jacobiean A_f of f with respect to u  and Q_f is an approximation of
754       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'
755       Q_g can be difficult, non-linear conjugate gradient method is applied to solve for p, so Q_g plays
756       in fact the role of a preconditioner.
757       """
758       def __init__(self,verbose=False,*args):
759           """
760           initializes the problem
761    
762           @param verbose: switches on the printing out some information
763           @type verbose: C{bool}
764           @note: this method may be overwritten by a particular saddle point problem
765           """
766           if not isinstance(verbose,bool):
767                raise TypeError("verbose needs to be of type bool.")
768           self.__verbose=verbose
769           self.relaxation=1.
770    
771       def trace(self,text):
772           """
773           prints text if verbose has been set
774    
775           @param text: a text message
776           @type text: C{str}
777           """
778           if self.__verbose: print "%s: %s"%(str(self),text)
779    
780       def solve_f(self,u,p,tol=1.e-8):
781           """
782           solves
783    
784           A_f du = f(u,p)
785    
786           with tolerance C{tol} and return du. A_f is Jacobiean of f with respect to u.
787    
788           @param u: current approximation of u
789           @type u: L{escript.Data}
790           @param p: current approximation of p
791           @type p: L{escript.Data}
792           @param tol: tolerance expected for du
793           @type tol: C{float}
794           @return: increment du
795           @rtype: L{escript.Data}
796           @note: this method has to be overwritten by a particular saddle point problem
797           """
798           pass
799    
800       def solve_g(self,u,tol=1.e-8):
801           """
802           solves
803    
804           Q_g dp = g(u)
805    
806           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.
807    
808           @param u: current approximation of u
809           @type u: L{escript.Data}
810           @param tol: tolerance expected for dp
811           @type tol: C{float}
812           @return: increment dp
813           @rtype: L{escript.Data}
814           @note: this method has to be overwritten by a particular saddle point problem
815           """
816           pass
817    
818       def inner(self,p0,p1):
819           """
820           inner product of p0 and p1 approximating p. Typically this returns integrate(p0*p1)
821           @return: inner product of p0 and p1
822           @rtype: C{float}
823           """
824           pass
825    
826       subiter_max=3
827       def solve(self,u0,p0,tolerance=1.e-6,tolerance_u=None,iter_max=100,accepted_reduction=0.995,relaxation=None):
828            """
829            runs the solver
830    
831            @param u0: initial guess for C{u}
832            @type u0: L{esys.escript.Data}
833            @param p0: initial guess for C{p}
834            @type p0: L{esys.escript.Data}
835            @param tolerance: tolerance for relative error in C{u} and C{p}
836            @type tolerance: positive C{float}
837            @param tolerance_u: tolerance for relative error in C{u} if different from C{tolerance}
838            @type tolerance_u: positive C{float}
839            @param iter_max: maximum number of iteration steps.
840            @type iter_max: C{int}
841            @param accepted_reduction: if the norm  g cannot be reduced by C{accepted_reduction} backtracking to adapt the
842                                       relaxation factor. If C{accepted_reduction=None} no backtracking is used.
843            @type accepted_reduction: positive C{float} or C{None}
844            @param relaxation: initial relaxation factor. If C{relaxation==None}, the last relaxation factor is used.
845            @type relaxation: C{float} or C{None}
846            """
847            tol=1.e-2
848            if tolerance_u==None: tolerance_u=tolerance
849            if not relaxation==None: self.relaxation=relaxation
850            if accepted_reduction ==None:
851                  angle_limit=0.
852            elif accepted_reduction>=1.:
853                  angle_limit=0.
854            else:
855                  angle_limit=util.sqrt(1-accepted_reduction**2)
856            self.iter=0
857            u=u0
858            p=p0
859            #
860            #   initialize things:
861            #
862            converged=False
863            #
864            #  start loop:
865            #
866            #  initial search direction is g
867            #
868            while not converged :
869                if self.iter>iter_max:
870                    raise ArithmeticError("no convergence after %s steps."%self.iter)
871                f_new=self.solve_f(u,p,tol)
872                norm_f_new = util.Lsup(f_new)
873                u_new=u-f_new
874                g_new=self.solve_g(u_new,tol)
875                self.iter+=1
876                norm_g_new = util.sqrt(self.inner(g_new,g_new))
877                if norm_f_new==0. and norm_g_new==0.: return u, p
878                if self.iter>1 and not accepted_reduction==None:
879                   #
880                   #   did we manage to reduce the norm of G? I
881                   #   if not we start a backtracking procedure
882                   #
883                   # print "new/old norm = ",norm_g_new, norm_g, norm_g_new/norm_g
884                   if norm_g_new > accepted_reduction * norm_g:
885                      sub_iter=0
886                      s=self.relaxation
887                      d=g
888                      g_last=g
889                      self.trace("    start substepping: f = %s, g = %s, relaxation = %s."%(norm_f_new, norm_g_new, s))
890                      while sub_iter < self.subiter_max and  norm_g_new > accepted_reduction * norm_g:
891                         dg= g_new-g_last
892                         norm_dg=abs(util.sqrt(self.inner(dg,dg))/self.relaxation)
893                         rad=self.inner(g_new,dg)/self.relaxation
894                         # print "   ",sub_iter,": rad, norm_dg:",abs(rad), norm_dg*norm_g_new * angle_limit
895                         # print "   ",sub_iter,": rad, norm_dg:",rad, norm_dg, norm_g_new, norm_g
896                         if abs(rad) < norm_dg*norm_g_new * angle_limit:
897                             if sub_iter>0: self.trace("    no further improvements expected from backtracking.")
898                             break
899                         r=self.relaxation
900                         self.relaxation= - rad/norm_dg**2
901                         s+=self.relaxation
902                         #####
903                         # a=g_new+self.relaxation*dg/r
904                         # print "predicted new norm = ",util.sqrt(self.inner(a,a)),util.sqrt(self.inner(g_new,g_new)), self.relaxation
905                         #####
906                         g_last=g_new
907                         p+=self.relaxation*d
908                         f_new=self.solve_f(u,p,tol)
909                         u_new=u-f_new
910                         g_new=self.solve_g(u_new,tol)
911                         self.iter+=1
912                         norm_f_new = util.Lsup(f_new)
913                         norm_g_new = util.sqrt(self.inner(g_new,g_new))
914                         # print "   ",sub_iter," new g norm",norm_g_new
915                         self.trace("    %s th sub-step: f = %s, g = %s, relaxation = %s."%(sub_iter, norm_f_new, norm_g_new, s))
916                         #
917                         #   can we expect reduction of g?
918                         #
919                         # u_last=u_new
920                         sub_iter+=1
921                      self.relaxation=s
922                #
923                #  check for convergence:
924                #
925                norm_u_new = util.Lsup(u_new)
926                p_new=p+self.relaxation*g_new
927                norm_p_new = util.sqrt(self.inner(p_new,p_new))
928                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))
929    
930                if self.iter>1:
931                   dg2=g_new-g
932                   df2=f_new-f
933                   norm_dg2=util.sqrt(self.inner(dg2,dg2))
934                   norm_df2=util.Lsup(df2)
935                   # print norm_g_new, norm_g, norm_dg, norm_p, tolerance
936                   tol_eq_g=tolerance*norm_dg2/(norm_g*abs(self.relaxation))*norm_p_new
937                   tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new
938                   if norm_g_new <= tol_eq_g and norm_f_new <= tol_eq_f:
939                       converged=True
940                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
941            self.trace("convergence after %s steps."%self.iter)
942            return u,p
943    #   def solve(self,u0,p0,tolerance=1.e-6,iter_max=10,self.relaxation=1.):
944    #      tol=1.e-2
945    #      iter=0
946    #      converged=False
947    #      u=u0*1.
948    #      p=p0*1.
949    #      while not converged and iter<iter_max:
950    #          du=self.solve_f(u,p,tol)
951    #          u-=du
952    #          norm_du=util.Lsup(du)
953    #          norm_u=util.Lsup(u)
954    #        
955    #          dp=self.relaxation*self.solve_g(u,tol)
956    #          p+=dp
957    #          norm_dp=util.sqrt(self.inner(dp,dp))
958    #          norm_p=util.sqrt(self.inner(p,p))
959    #          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)
960    #          iter+=1
961    #
962    #          converged = (norm_du <= tolerance*norm_u) and  (norm_dp <= tolerance*norm_p)
963    #      if converged:
964    #          print "convergence after %s steps."%iter
965    #      else:
966    #          raise ArithmeticError("no convergence after %s steps."%iter)
967    #
968    #      return u,p
969              
970    def MaskFromBoundaryTag(function_space,*tags):
971       """
972       create a mask on the given function space which one for samples
973       that touch the boundary tagged by tags.
974    
975       usage: m=MaskFromBoundaryTag(Solution(domain),"left", "right")
976    
977       @param function_space: a given function space
978       @type function_space: L{escript.FunctionSpace}
979       @param tags: boundray tags
980       @type tags: C{str}
981       @return: a mask which marks samples used by C{function_space} that are touching the
982                boundary tagged by any of the given tags.
983       @rtype: L{escript.Data} of rank 0
984       """
985       pde=linearPDEs.LinearPDE(function_space.getDomain(),numEquations=1, numSolutions=1)
986       d=escript.Scalar(0.,escript.FunctionOnBoundary(function_space.getDomain()))
987       for t in tags: d.setTaggedValue(t,1.)
988       pde.setValue(y=d)
989       out=util.whereNonZero(pde.getRightHandSide())
990       if out.getFunctionSpace() == function_space:
991          return out
992       else:
993          return util.whereNonZero(util.interpolate(out,function_space))
994    
995    

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

  ViewVC Help
Powered by ViewVC 1.1.26