/[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 614 by elspeth, Wed Mar 22 01:37:07 2006 UTC revision 1519 by artak, Tue Apr 22 03:45:36 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
27    @var __copyright__: copyrights
28    @var __license__: licence agreement
29    @var __url__: url entry point on documentation
30    @var __version__: version
31    @var __date__: date of the version
32  """  """
33    
34    __author__="Lutz Gross, l.gross@uq.edu.au"
35  __copyright__="""  Copyright (c) 2006 by ACcESS MNRF  __copyright__="""  Copyright (c) 2006 by ACcESS MNRF
36                      http://www.access.edu.au                      http://www.access.edu.au
37                  Primary Business: Queensland, Australia"""                  Primary Business: Queensland, Australia"""
38  __license__="""Licensed under the Open Software License version 3.0  __license__="""Licensed under the Open Software License version 3.0
39               http://www.opensource.org/licenses/osl-3.0.php"""               http://www.opensource.org/licenses/osl-3.0.php"""
40    __url__="http://www.iservo.edu.au/esys"
41    __version__="$Revision$"
42    __date__="$Date$"
43    
44    
45  import escript  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    """    """
58    a simple mechanism to manage time dependend values.    a simple mechanism to manage time dependend values.
59    
60    typical usage is:    typical usage is::
61    
62    dt=0.1 # time increment       dt=0.1 # time increment
63    tm=TimeIntegrationManager(inital_value,p=1)       tm=TimeIntegrationManager(inital_value,p=1)
64    while t<1.       while t<1.
65        v_guess=tm.extrapolate(dt) # extrapolate to t+dt           v_guess=tm.extrapolate(dt) # extrapolate to t+dt
66        v=...           v=...
67        tm.checkin(dt,v)           tm.checkin(dt,v)
68        t+=dt           t+=dt
69    
70    @remark: currently only p=1 is supported.    @note: currently only p=1 is supported.
71    """    """
72    def __init__(self,*inital_values,**kwargs):    def __init__(self,*inital_values,**kwargs):
73       """       """
# Line 119  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 129  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 180  class NoPDE: Line 211  class NoPDE:
211    
212       The constraint is overwriting any other condition.       The constraint is overwriting any other condition.
213    
214       @remark: This class is similar to the L{LinearPDE} class with A=B=C=X=0 but has the intention       @note: This class is similar to the L{linearPDEs.LinearPDE} class with A=B=C=X=0 but has the intention
215                that all input parameter are given in L{Solution} or L{ReducedSolution}. The whole              that all input parameter are given in L{Solution} or L{ReducedSolution}. The whole
216                thing is a bit strange and I blame Robert.Woodcock@csiro.au for this.              thing is a bit strange and I blame Robert.Woodcock@csiro.au for this.
217       """       """
218       def __init__(self,domain,D=None,Y=None,q=None,r=None):       def __init__(self,domain,D=None,Y=None,q=None,r=None):
219           """           """
# Line 191  class NoPDE: Line 222  class NoPDE:
222           @param domain: domain of the PDE.           @param domain: domain of the PDE.
223           @type domain: L{Domain}           @type domain: L{Domain}
224           @param D: coefficient of the solution.           @param D: coefficient of the solution.
225           @type D: C{float}, C{int}, L{NumArray}, L{Data}           @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}
226           @param Y: right hand side           @param Y: right hand side
227           @type Y: C{float}, C{int}, L{NumArray}, L{Data}           @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}
228           @param q: location of constraints           @param q: location of constraints
229           @type q: C{float}, C{int}, L{NumArray}, L{Data}           @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}
230           @param r: value of solution at locations of constraints           @param r: value of solution at locations of constraints
231           @type r: C{float}, C{int}, L{NumArray}, L{Data}           @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}
232           """           """
233           self.__domain=domain           self.__domain=domain
234           self.__D=D           self.__D=D
# Line 225  class NoPDE: Line 256  class NoPDE:
256           assigns values to the parameters.           assigns values to the parameters.
257    
258           @param D: coefficient of the solution.           @param D: coefficient of the solution.
259           @type D: C{float}, C{int}, L{NumArray}, L{Data}           @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}
260           @param Y: right hand side           @param Y: right hand side
261           @type Y: C{float}, C{int}, L{NumArray}, L{Data}           @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}
262           @param q: location of constraints           @param q: location of constraints
263           @type q: C{float}, C{int}, L{NumArray}, L{Data}           @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}
264           @param r: value of solution at locations of constraints           @param r: value of solution at locations of constraints
265           @type r: C{float}, C{int}, L{NumArray}, L{Data}           @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}
266           """           """
267           if not D==None:           if not D==None:
268              self.__D=D              self.__D=D
# Line 284  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 303  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 329  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
394             else:             else:
395               out=data.interpolate(self.getFunctionSpace()).convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])               dat=data.interpolate(self.getFunctionSpace())
396             if data.getRank()==0:             id=self.getId()
397                return out[0]             r=data.getRank()
398               if isinstance(id,list):
399                   out=[]
400                   for i in id:
401                      o=data.getValueOfGlobalDataPoint(*i)
402                      if data.getRank()==0:
403                         out.append(o[0])
404                      else:
405                         out.append(o)
406                   return out
407             else:             else:
408                return out               out=data.getValueOfGlobalDataPoint(*id)
409                 if data.getRank()==0:
410                    return out[0]
411                 else:
412                    return out
413          else:          else:
414             return data             return data
415    
416  # vim: expandtab shiftwidth=4:  class SolverSchemeException(Exception):
417       """
418       exceptions thrown by solvers
419       """
420       pass
421    
422    class IndefinitePreconditioner(SolverSchemeException):
423       """
424       the preconditioner is not positive definite.
425       """
426       pass
427    class MaxIterReached(SolverSchemeException):
428       """
429       maxium number of iteration steps is reached.
430       """
431       pass
432    class IterationBreakDown(SolverSchemeException):
433       """
434       iteration scheme econouters an incurable breakdown.
435       """
436       pass
437    class NegativeNorm(SolverSchemeException):
438       """
439       a norm calculation returns a negative norm.
440       """
441       pass
442    
443    class IterationHistory(object):
444       """
445       The IterationHistory class is used to define a stopping criterium. It keeps track of the
446       residual norms. The stoppingcriterium indicates termination if the residual norm has been reduced by
447       a given tolerance.
448       """
449       def __init__(self,tolerance=math.sqrt(util.EPSILON),verbose=False):
450          """
451          Initialization
452    
453          @param tolerance: tolerance
454          @type tolerance: positive C{float}
455          @param verbose: switches on the printing out some information
456          @type verbose: C{bool}
457          """
458          if not tolerance>0.:
459              raise ValueError,"tolerance needs to be positive."
460          self.tolerance=tolerance
461          self.verbose=verbose
462          self.history=[]
463       def stoppingcriterium(self,norm_r,r,x):
464           """
465           returns True if the C{norm_r} is C{tolerance}*C{norm_r[0]} where C{norm_r[0]}  is the residual norm at the first call.
466    
467          
468           @param norm_r: current residual norm
469           @type norm_r: non-negative C{float}
470           @param r: current residual (not used)
471           @param x: current solution approximation (not used)
472           @return: C{True} is the stopping criterium is fullfilled. Otherwise C{False} is returned.
473           @rtype: C{bool}
474    
475           """
476           self.history.append(norm_r)
477           if self.verbose: print "iter: #s:  inner(rhat,r) = #e"#(len(self.history)-1, self.history[-1])
478           return self.history[-1]<=self.tolerance * self.history[0]
479    
480       def stoppingcriterium2(self,norm_r,norm_b,solver="GMRES",TOL=None):
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           if TOL==None:
494              TOL=self.tolerance
495           self.history.append(norm_r)
496           if self.verbose: print "iter: #s:  norm(r) = #e"#(len(self.history)-1, self.history[-1])
497           return self.history[-1]<=TOL * norm_b
498    
499    def PCG(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
500       """
501       Solver for
502    
503       M{Ax=b}
504    
505       with a symmetric and positive definite operator A (more details required!).
506       It uses the conjugate gradient method with preconditioner M providing an approximation of A.
507    
508       The iteration is terminated if the C{stoppingcriterium} function return C{True}.
509    
510       For details on the preconditioned conjugate gradient method see the book:
511    
512       Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
513       T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
514       C. Romine, and H. van der Vorst.
515    
516       @param b: the right hand side of the liner system. C{b} is altered.
517       @type b: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
518       @param Aprod: returns the value Ax
519       @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}.
520       @param Msolve: solves Mx=r
521       @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
522    type like argument C{x}.
523       @param bilinearform: inner product C{<x,r>}
524       @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}.
525       @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.
526       @type stoppingcriterium: function that returns C{True} or C{False}
527       @param x: an initial guess for the solution. If no C{x} is given 0*b is used.
528       @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
529       @param iter_max: maximum number of iteration steps.
530       @type iter_max: C{int}
531       @return: the solution approximation and the corresponding residual
532       @rtype: C{tuple}
533       @warning: C{b} and C{x} are altered.
534       """
535       iter=0
536       if x==None:
537          x=0*b
538       else:
539          b += (-1)*Aprod(x)
540       r=b
541       rhat=Msolve(r)
542       d = rhat
543       rhat_dot_r = bilinearform(rhat, r)
544       if rhat_dot_r<0: raise NegativeNorm,"negative norm."
545    
546       while not stoppingcriterium(math.sqrt(rhat_dot_r),r,x):
547           iter+=1
548           if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
549    
550           q=Aprod(d)
551           alpha = rhat_dot_r / bilinearform(d, q)
552           x += alpha * d
553           r += (-alpha) * q
554    
555           rhat=Msolve(r)
556           rhat_dot_r_new = bilinearform(rhat, r)
557           beta = rhat_dot_r_new / rhat_dot_r
558           rhat+=beta * d
559           d=rhat
560    
561           rhat_dot_r = rhat_dot_r_new
562           if rhat_dot_r<0: raise NegativeNorm,"negative norm."
563    
564       return x,r
565    
566    
567    ############################
568    # Added by Artak
569    #################################3
570    
571    #Apply a sequence of k Givens rotations, used within gmres codes
572    # vrot=givapp(c, s, vin, k)
573    def givapp(c,s,vin):
574        vrot=vin # warning: vin is altered!!!!
575        if isinstance(c,float):
576            vrot=[c*vrot[0]-s*vrot[1],s*vrot[0]+c*vrot[1]]
577        else:
578            for i in range(len(c)):
579                w1=c[i]*vrot[i]-s[i]*vrot[i+1]
580            w2=s[i]*vrot[i]+c[i]*vrot[i+1]
581                vrot[i:i+2]=w1,w2
582        return vrot
583    
584    def GMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):
585       m=iter_restart
586       iter=0
587       while True:
588          if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached"%iter_max
589          x,stopped=GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=x, iter_max=iter_max-iter, iter_restart=m)
590          if stopped: break
591          iter+=iter_restart    
592       return x
593    
594    def GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):
595       iter=0
596       r=Msolve(b)
597       r_dot_r = bilinearform(r, r)
598       if r_dot_r<0: raise NegativeNorm,"negative norm."
599       norm_b=math.sqrt(r_dot_r)
600    
601       if x==None:
602          x=0
603       else:
604          r=Msolve(b-Aprod(x))
605          r_dot_r = bilinearform(r, r)
606          if r_dot_r<0: raise NegativeNorm,"negative norm."
607      
608       h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)
609       c=numarray.zeros(iter_restart,numarray.Float64)
610       s=numarray.zeros(iter_restart,numarray.Float64)
611       g=numarray.zeros(iter_restart,numarray.Float64)
612       v=[]
613    
614       rho=math.sqrt(r_dot_r)
615      
616       v.append(r/rho)
617       g[0]=rho
618       while not (stoppingcriterium(rho,norm_b) or iter==iter_restart-1):
619    
620        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
621    
622        
623        p=Msolve(Aprod(v[iter]))
624    
625        v.append(p)
626    
627        v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
628    
629    # Modified Gram-Schmidt
630        for j in range(iter+1):
631          h[j][iter]=bilinearform(v[j],v[iter+1])  
632          v[iter+1]+=(-1.)*h[j][iter]*v[j]
633          
634        h[iter+1][iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
635        v_norm2=h[iter+1][iter]
636    
637    
638    # Reorthogonalize if needed
639        if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
640         for j in range(iter+1):
641            hr=bilinearform(v[j],v[iter+1])
642                h[j][iter]=h[j][iter]+hr #vhat
643                v[iter+1] +=(-1.)*hr*v[j]
644    
645         v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
646         h[iter+1][iter]=v_norm2
647    
648    #   watch out for happy breakdown
649            if v_norm2 != 0:
650             v[iter+1]=v[iter+1]/h[iter+1][iter]
651    
652    #   Form and store the information for the new Givens rotation
653        if iter > 0 :
654            hhat=[]
655            for i in range(iter+1) : hhat.append(h[i][iter])
656            hhat=givapp(c[0:iter],s[0:iter],hhat);
657                for i in range(iter+1) : h[i][iter]=hhat[i]
658    
659        mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])
660        if mu!=0 :
661            c[iter]=h[iter][iter]/mu
662            s[iter]=-h[iter+1][iter]/mu
663            h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter]
664            h[iter+1][iter]=0.0
665            g[iter:iter+2]=givapp(c[iter],s[iter],g[iter:iter+2])
666    
667    # Update the residual norm
668            rho=abs(g[iter+1])
669        iter+=1
670    
671    # At this point either iter > iter_max or rho < tol.
672    # It's time to compute x and leave.        
673    
674       if iter > 0 :
675         y=numarray.zeros(iter,numarray.Float64)    
676         y[iter-1] = g[iter-1] / h[iter-1][iter-1]
677         if iter > 1 :  
678            i=iter-2  
679            while i>=0 :
680              y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]
681              i=i-1
682         xhat=v[iter-1]*y[iter-1]
683         for i in range(iter-1):
684        xhat += v[i]*y[i]
685       else : xhat=v[0]
686        
687       x += xhat
688       if iter<iter_restart-1:
689          stopped=True
690       else:
691          stopped=False
692    
693       return x,stopped
694    
695    ######################################################
696    def dirder(x, w, bilinearform, Aprod, Msolve, f0, b ):
697    ######################################################
698    
699    # DIRDER estimates the directional derivative of a function F.
700    
701    
702    # Hardwired difference increment.
703    #
704      epsnew = 1.0e-07
705    #
706    #  Scale the step.
707    #
708      norm_w=math.sqrt(bilinearform(w,w))
709      if norm_w== 0.0:
710        return x/x
711    
712      epsnew = epsnew / norm_w
713    
714      if norm_w > 0.0:
715        epsnew = epsnew * math.sqrt(bilinearform(x,x))
716    #
717    #  DEL and F1 could share the same space if storage
718    #  is more important than clarity.
719    #
720    
721      DEL = x + epsnew * w
722      f1 = -Msolve(Aprod(DEL))
723      z = ( f1 - f0 ) / epsnew
724      return z
725    
726    ######################################################
727    def FDGMRES(f0, Aprod, Msolve, bilinearform, stoppingcriterium, xc=None, x=None, iter_max=100, iter_restart=20,TOL=None):
728    ######################################################
729       b=-f0
730       b_dot_b = bilinearform(b, b)
731       if b_dot_b<0: raise NegativeNorm,"negative norm."
732       norm_b=math.sqrt(b_dot_b)
733    
734       r=b
735    
736       if x==None:
737          x=0
738       else:
739          r=-dirder(xc,x,bilinearform,Aprod,Msolve,f0,b)-f0  
740          
741       r_dot_r = bilinearform(r, r)
742       if r_dot_r<0: raise NegativeNorm,"negative norm."
743      
744       h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)
745       c=numarray.zeros(iter_restart,numarray.Float64)
746       s=numarray.zeros(iter_restart,numarray.Float64)
747       g=numarray.zeros(iter_restart,numarray.Float64)
748       v=[]
749    
750       rho=math.sqrt(r_dot_r)
751      
752       v.append(r/rho)
753       g[0]=rho
754       iter=0
755    
756       while not (stoppingcriterium(rho,norm_b,solver="FDGMRES",TOL=TOL) or iter==iter_restart-1):
757    
758        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
759    
760        
761            p=dirder(xc, v[iter], bilinearform,Aprod,Msolve,f0,b)
762    
763        v.append(p)
764    
765        v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
766    
767    # Modified Gram-Schmidt
768        for j in range(iter+1):
769          h[j][iter]=bilinearform(v[j],v[iter+1])  
770          v[iter+1]+=(-1.)*h[j][iter]*v[j]
771          
772        h[iter+1][iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
773        v_norm2=h[iter+1][iter]
774    
775    
776    # Reorthogonalize if needed
777        if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
778         for j in range(iter+1):
779            hr=bilinearform(v[j],v[iter+1])
780                h[j][iter]=h[j][iter]+hr #vhat
781                v[iter+1] +=(-1.)*hr*v[j]
782    
783         v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
784         h[iter+1][iter]=v_norm2
785    
786    #   watch out for happy breakdown
787            if v_norm2 != 0:
788             v[iter+1]=v[iter+1]/h[iter+1][iter]
789    
790    #   Form and store the information for the new Givens rotation
791        if iter > 0 :
792            hhat=[]
793            for i in range(iter+1) : hhat.append(h[i][iter])
794            hhat=givapp(c[0:iter],s[0:iter],hhat);
795                for i in range(iter+1) : h[i][iter]=hhat[i]
796    
797        mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])
798        if mu!=0 :
799            c[iter]=h[iter][iter]/mu
800            s[iter]=-h[iter+1][iter]/mu
801            h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter]
802            h[iter+1][iter]=0.0
803            g[iter:iter+2]=givapp(c[iter],s[iter],g[iter:iter+2])
804    
805    # Update the residual norm
806            rho=abs(g[iter+1])
807        iter+=1
808    
809       if iter > 0 :
810         y=numarray.zeros(iter,numarray.Float64)    
811         y[iter-1] = g[iter-1] / h[iter-1][iter-1]
812         if iter > 1 :  
813            i=iter-2  
814            while i>=0 :
815              y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]
816              i=i-1
817         xhat=v[iter-1]*y[iter-1]
818         for i in range(iter-1):
819        xhat += v[i]*y[i]
820       else : xhat=v[0]
821        
822       x += xhat
823       if iter<iter_restart-1:
824          stopped=True
825       else:
826          stopped=False
827    
828       return x,stopped
829    
830    def MINRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
831    
832        #
833        #  minres solves the system of linear equations Ax = b
834        #  where A is a symmetric matrix (possibly indefinite or singular)
835        #  and b is a given vector.
836        #  
837        #  "A" may be a dense or sparse matrix (preferably sparse!)
838        #  or the name of a function such that
839        #               y = A(x)
840        #  returns the product y = Ax for any given vector x.
841        #
842        #  "M" defines a positive-definite preconditioner M = C C'.
843        #  "M" may be a dense or sparse matrix (preferably sparse!)
844        #  or the name of a function such that
845        #  solves the system My = x for any given vector x.
846        #
847        #
848        
849        #------------------------------------------------------------------
850        # Set up y and v for the first Lanczos vector v1.
851        # y  =  beta1 P' v1,  where  P = C**(-1).
852        # v is really P' v1.
853        #------------------------------------------------------------------
854        if x==None:
855          x=0*b
856        else:
857          b += (-1)*Aprod(x)
858    
859        r1    = b
860        y = Msolve(b)
861        beta1 = bilinearform(b,y)
862    
863        if beta1< 0: raise NegativeNorm,"negative norm."
864    
865        #  If b = 0 exactly, stop with x = 0.
866        if beta1==0: return x*0.
867    
868        if beta1> 0:
869          beta1  = math.sqrt(beta1)      
870    
871        #------------------------------------------------------------------
872        # Initialize quantities.
873        # ------------------------------------------------------------------
874        iter   = 0
875        Anorm = 0
876        ynorm = 0
877        oldb   = 0
878        beta   = beta1
879        dbar   = 0
880        epsln  = 0
881        phibar = beta1
882        rhs1   = beta1
883        rhs2   = 0
884        rnorm  = phibar
885        tnorm2 = 0
886        ynorm2 = 0
887        cs     = -1
888        sn     = 0
889        w      = b*0.
890        w2     = b*0.
891        r2     = r1
892        eps    = 0.0001
893    
894        #---------------------------------------------------------------------
895        # Main iteration loop.
896        # --------------------------------------------------------------------
897        while not stoppingcriterium(rnorm,Anorm*ynorm,'MINRES'):    #  checks ||r|| < (||A|| ||x||) * TOL
898    
899        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
900            iter    = iter  +  1
901    
902            #-----------------------------------------------------------------
903            # Obtain quantities for the next Lanczos vector vk+1, k = 1, 2,...
904            # The general iteration is similar to the case k = 1 with v0 = 0:
905            #
906            #   p1      = Operator * v1  -  beta1 * v0,
907            #   alpha1  = v1'p1,
908            #   q2      = p2  -  alpha1 * v1,
909            #   beta2^2 = q2'q2,
910            #   v2      = (1/beta2) q2.
911            #
912            # Again, y = betak P vk,  where  P = C**(-1).
913            #-----------------------------------------------------------------
914            s = 1/beta                 # Normalize previous vector (in y).
915            v = s*y                    # v = vk if P = I
916        
917            y      = Aprod(v)
918        
919            if iter >= 2:
920              y = y - (beta/oldb)*r1
921    
922            alfa   = bilinearform(v,y)              # alphak
923            y      = (- alfa/beta)*r2 + y
924            r1     = r2
925            r2     = y
926            y = Msolve(r2)
927            oldb   = beta                           # oldb = betak
928            beta   = bilinearform(r2,y)             # beta = betak+1^2
929            if beta < 0: raise NegativeNorm,"negative norm."
930    
931            beta   = math.sqrt( beta )
932            tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta
933            
934            if iter==1:                 # Initialize a few things.
935              gmax   = abs( alfa )      # alpha1
936              gmin   = gmax             # alpha1
937    
938            # Apply previous rotation Qk-1 to get
939            #   [deltak epslnk+1] = [cs  sn][dbark    0   ]
940            #   [gbar k dbar k+1]   [sn -cs][alfak betak+1].
941        
942            oldeps = epsln
943            delta  = cs * dbar  +  sn * alfa  # delta1 = 0         deltak
944            gbar   = sn * dbar  -  cs * alfa  # gbar 1 = alfa1     gbar k
945            epsln  =               sn * beta  # epsln2 = 0         epslnk+1
946            dbar   =            -  cs * beta  # dbar 2 = beta2     dbar k+1
947    
948            # Compute the next plane rotation Qk
949    
950            gamma  = math.sqrt(gbar*gbar+beta*beta)  # gammak
951            gamma  = max(gamma,eps)
952            cs     = gbar / gamma             # ck
953            sn     = beta / gamma             # sk
954            phi    = cs * phibar              # phik
955            phibar = sn * phibar              # phibark+1
956    
957            # Update  x.
958    
959            denom = 1/gamma
960            w1    = w2
961            w2    = w
962            w     = (v - oldeps*w1 - delta*w2) * denom
963            x     = x  +  phi*w
964    
965            # Go round again.
966    
967            gmax   = max(gmax,gamma)
968            gmin   = min(gmin,gamma)
969            z      = rhs1 / gamma
970            ynorm2 = z*z  +  ynorm2
971            rhs1   = rhs2 -  delta*z
972            rhs2   =      -  epsln*z
973    
974            # Estimate various norms and test for convergence.
975    
976            Anorm  = math.sqrt( tnorm2 )
977            ynorm  = math.sqrt( ynorm2 )
978    
979            rnorm  = phibar
980    
981        return x
982        
983    def NewtonGMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium,x=None, iter_max=100,iter_restart=20):
984    
985        gamma=.9
986        lmaxit=40
987        etamax=.5
988    
989        n = 1 #len(x)
990        iter=0
991        
992        # evaluate f at the initial iterate
993        # compute the stop tolerance
994        #
995        r=b
996        if x==None:
997          x=0*b
998        else:
999          r += (-1)*Aprod(x)
1000    
1001        f0=-Msolve(r)
1002        fnrm=math.sqrt(bilinearform(f0,f0))/math.sqrt(n)
1003        fnrmo=1
1004        atol=1.e-2
1005        rtol=1.e-4
1006        stop_tol=atol + rtol*fnrm
1007        #
1008        # main iteration loop
1009        #
1010        while not stoppingcriterium(fnrm,stop_tol,'NewtonGMRES',TOL=1.):
1011    
1012                if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
1013            #
1014            # keep track of the ratio (rat = fnrm/frnmo)
1015            # of successive residual norms and
1016            # the iteration counter (iter)
1017            #
1018            #rat=fnrm/fnrmo
1019            fnrmo=fnrm
1020            iter+=1
1021            #
1022                # compute the step using a GMRES(m) routine especially designed for this purpose
1023            #
1024                initer=0
1025                while True:
1026                   xc,stopped=FDGMRES(f0, Aprod, Msolve, bilinearform, stoppingcriterium, xc=x, iter_max=lmaxit-initer, iter_restart=iter_restart, TOL=etamax)
1027                   if stopped: break
1028                   initer+=iter_restart
1029            xold=x
1030            x+=xc
1031            f0=-Msolve(Aprod(x))
1032            fnrm=math.sqrt(bilinearform(f0,f0))/math.sqrt(n)
1033            rat=fnrm/fnrmo
1034    
1035    
1036        #   adjust eta
1037        #
1038            if etamax > 0:
1039                etaold=etamax
1040                etanew=gamma*rat*rat
1041                if gamma*etaold*etaold > .1 :
1042                    etanew=max(etanew,gamma*etaold*etaold)
1043                
1044                etamax=min(etanew,etamax)
1045                etamax=max(etamax,.5*stop_tol/fnrm)
1046    
1047        return x
1048    
1049    def TFQMR(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
1050    
1051    # TFQMR solver for linear systems
1052    #
1053    #
1054    # initialization
1055    #
1056      errtol = math.sqrt(bilinearform(b,b))
1057      norm_b=errtol
1058      kmax  = iter_max
1059      error = []
1060    
1061      if math.sqrt(bilinearform(x,x)) != 0.0:
1062        r = b - Aprod(x)
1063      else:
1064        r = b
1065    
1066      r=Msolve(r)
1067    
1068      u1=0
1069      u2=0
1070      y1=0
1071      y2=0
1072    
1073      w = r
1074      y1 = r
1075      iter = 0
1076      d = 0
1077      
1078      v = Msolve(Aprod(y1))
1079      u1 = v
1080      
1081      theta = 0.0;
1082      eta = 0.0;
1083      tau = math.sqrt(bilinearform(r,r))
1084      error = [ error, tau ]
1085      rho = tau * tau
1086      m=1
1087    #
1088    #  TFQMR iteration
1089    #
1090    #  while ( iter < kmax-1 ):
1091      
1092      while not stoppingcriterium(tau*math.sqrt ( m + 1 ),norm_b,'TFQMR'):
1093        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
1094    
1095        sigma = bilinearform(r,v)
1096    
1097        if ( sigma == 0.0 ):
1098          raise 'TFQMR breakdown, sigma=0'
1099        
1100    
1101        alpha = rho / sigma
1102    
1103        for j in range(2):
1104    #
1105    #   Compute y2 and u2 only if you have to
1106    #
1107          if ( j == 1 ):
1108            y2 = y1 - alpha * v
1109            u2 = Msolve(Aprod(y2))
1110    
1111          m = 2 * (iter+1) - 2 + (j+1)
1112          if j==0:
1113             w = w - alpha * u1
1114             d = y1 + ( theta * theta * eta / alpha ) * d
1115          if j==1:
1116             w = w - alpha * u2
1117             d = y2 + ( theta * theta * eta / alpha ) * d
1118    
1119          theta = math.sqrt(bilinearform(w,w))/ tau
1120          c = 1.0 / math.sqrt ( 1.0 + theta * theta )
1121          tau = tau * theta * c
1122          eta = c * c * alpha
1123          x = x + eta * d
1124    #
1125    #  Try to terminate the iteration at each pass through the loop
1126    #
1127         # if ( tau * math.sqrt ( m + 1 ) <= errtol ):
1128         #   error = [ error, tau ]
1129         #   total_iters = iter
1130         #   break
1131          
1132    
1133        if ( rho == 0.0 ):
1134          raise 'TFQMR breakdown, rho=0'
1135        
1136    
1137        rhon = bilinearform(r,w)
1138        beta = rhon / rho;
1139        rho = rhon;
1140        y1 = w + beta * y2;
1141        u1 = Msolve(Aprod(y1))
1142        v = u1 + beta * ( u2 + beta * v )
1143        error = [ error, tau ]
1144        total_iters = iter
1145        
1146        iter = iter + 1
1147    
1148      return x
1149    
1150    
1151    #############################################
1152    
1153    class ArithmeticTuple(object):
1154       """
1155       tuple supporting inplace update x+=y and scaling x=a*y where x,y is an ArithmeticTuple and a is a float.
1156    
1157       example of usage:
1158    
1159       from esys.escript import Data
1160       from numarray import array
1161       a=Data(...)
1162       b=array([1.,4.])
1163       x=ArithmeticTuple(a,b)
1164       y=5.*x
1165    
1166       """
1167       def __init__(self,*args):
1168           """
1169           initialize object with elements args.
1170    
1171           @param args: tuple of object that support implace add (x+=y) and scaling (x=a*y)
1172           """
1173           self.__items=list(args)
1174    
1175       def __len__(self):
1176           """
1177           number of items
1178    
1179           @return: number of items
1180           @rtype: C{int}
1181           """
1182           return len(self.__items)
1183    
1184       def __getitem__(self,index):
1185           """
1186           get an item
1187    
1188           @param index: item to be returned
1189           @type index: C{int}
1190           @return: item with index C{index}
1191           """
1192           return self.__items.__getitem__(index)
1193    
1194       def __mul__(self,other):
1195           """
1196           scaling from the right
1197    
1198           @param other: scaling factor
1199           @type other: C{float}
1200           @return: itemwise self*other
1201           @rtype: L{ArithmeticTuple}
1202           """
1203           out=[]
1204           for i in range(len(self)):
1205               out.append(self[i]*other)
1206           return ArithmeticTuple(*tuple(out))
1207    
1208       def __rmul__(self,other):
1209           """
1210           scaling from the left
1211    
1212           @param other: scaling factor
1213           @type other: C{float}
1214           @return: itemwise other*self
1215           @rtype: L{ArithmeticTuple}
1216           """
1217           out=[]
1218           for i in range(len(self)):
1219               out.append(other*self[i])
1220           return ArithmeticTuple(*tuple(out))
1221    
1222    #########################
1223    # Added by Artak
1224    #########################
1225       def __div__(self,other):
1226           """
1227           dividing from the right
1228    
1229           @param other: scaling factor
1230           @type other: C{float}
1231           @return: itemwise self/other
1232           @rtype: L{ArithmeticTuple}
1233           """
1234           out=[]
1235           for i in range(len(self)):
1236               out.append(self[i]/other)
1237           return ArithmeticTuple(*tuple(out))
1238    
1239       def __rdiv__(self,other):
1240           """
1241           dividing from the left
1242    
1243           @param other: scaling factor
1244           @type other: C{float}
1245           @return: itemwise other/self
1246           @rtype: L{ArithmeticTuple}
1247           """
1248           out=[]
1249           for i in range(len(self)):
1250               out.append(other/self[i])
1251           return ArithmeticTuple(*tuple(out))
1252      
1253    ##########################################33
1254    
1255       def __iadd__(self,other):
1256           """
1257           in-place add of other to self
1258    
1259           @param other: increment
1260           @type other: C{ArithmeticTuple}
1261           """
1262           if len(self) != len(other):
1263               raise ValueError,"tuple length must match."
1264           for i in range(len(self)):
1265               self.__items[i]+=other[i]
1266           return self
1267    
1268    class HomogeneousSaddlePointProblem(object):
1269          """
1270          This provides a framwork for solving homogeneous saddle point problem of the form
1271    
1272                 Av+B^*p=f
1273                 Bv    =0
1274    
1275          for the unknowns v and p and given operators A and B and given right hand side f.
1276          B^* is the adjoint operator of B is the given inner product.
1277    
1278          """
1279          def __init__(self,**kwargs):
1280            self.setTolerance()
1281            self.setToleranceReductionFactor()
1282    
1283          def initialize(self):
1284            """
1285            initialize the problem (overwrite)
1286            """
1287            pass
1288          def B(self,v):
1289             """
1290             returns Bv (overwrite)
1291             @rtype: equal to the type of p
1292    
1293             @note: boundary conditions on p should be zero!
1294             """
1295             pass
1296    
1297          def inner(self,p0,p1):
1298             """
1299             returns inner product of two element p0 and p1  (overwrite)
1300            
1301             @type p0: equal to the type of p
1302             @type p1: equal to the type of p
1303             @rtype: C{float}
1304    
1305             @rtype: equal to the type of p
1306             """
1307             pass
1308    
1309          def solve_A(self,u,p):
1310             """
1311             solves Av=f-Au-B^*p with accuracy self.getReducedTolerance() (overwrite)
1312    
1313             @rtype: equal to the type of v
1314             @note: boundary conditions on v should be zero!
1315             """
1316             pass
1317    
1318          def solve_prec(self,p):
1319             """
1320             provides a preconditioner for BA^{-1}B^* with accuracy self.getReducedTolerance() (overwrite)
1321    
1322             @rtype: equal to the type of p
1323             """
1324             pass
1325    
1326          def stoppingcriterium(self,Bv,v,p):
1327             """
1328             returns a True if iteration is terminated. (overwrite)
1329    
1330             @rtype: C{bool}
1331             """
1332             pass
1333                
1334          def __inner(self,p,r):
1335             return self.inner(p,r[1])
1336    
1337          def __inner_p(self,p1,p2):
1338             return self.inner(p1,p2)
1339    
1340          def __stoppingcriterium(self,norm_r,r,p):
1341              return self.stoppingcriterium(r[1],r[0],p)
1342    
1343          def __stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None):
1344              return self.stoppingcriterium2(norm_r,norm_b,solver,TOL)
1345    
1346          def setTolerance(self,tolerance=1.e-8):
1347                  self.__tol=tolerance
1348          def getTolerance(self):
1349                  return self.__tol
1350          def setToleranceReductionFactor(self,reduction=0.01):
1351                  self.__reduction=reduction
1352          def getSubProblemTolerance(self):
1353                  return self.__reduction*self.getTolerance()
1354    
1355          def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='PCG',iter_restart=20):
1356                  """
1357                  solves the saddle point problem using initial guesses v and p.
1358    
1359                  @param max_iter: maximum number of iteration steps.
1360                  """
1361                  self.verbose=verbose
1362                  self.show_details=show_details and self.verbose
1363    
1364                  # assume p is known: then v=A^-1(f-B^*p)
1365                  # which leads to BA^-1B^*p = BA^-1f  
1366    
1367              # Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask)      
1368    
1369              
1370              self.__z=v+self.solve_A(v,p*0)
1371    
1372                  Bz=self.B(self.__z)
1373                  #
1374              #   solve BA^-1B^*p = Bz
1375                  #
1376                  #   note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv
1377                  #
1378                  #   with                    Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask)
1379                  #                           A(v-z)=Az-B^*p-Az = f -Az - B^*p (v-z=0 on fixed_u_mask)
1380                  #
1381                  self.iter=0
1382              if solver=='GMRES':      
1383                    if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter
1384                    p=GMRES(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.,iter_restart=iter_restart)
1385                    # solve Au=f-B^*p
1386                    #       A(u-v)=f-B^*p-Av
1387                    #       u=v+(u-v)
1388            u=v+self.solve_A(v,p)
1389    
1390              if solver=='NewtonGMRES':    
1391                    if self.verbose: print "enter NewtonGMRES method (iter_max=%s)"%max_iter
1392                    p=NewtonGMRES(Bz,self.__Aprod_Newton,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)
1393                    # solve Au=f-B^*p
1394                    #       A(u-v)=f-B^*p-Av
1395                    #       u=v+(u-v)
1396            u=v+self.solve_A(v,p)
1397                    
1398    
1399              if solver=='TFQMR':      
1400                    if self.verbose: print "enter TFQMR method (iter_max=%s)"%max_iter
1401                    p=TFQMR(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)
1402                    # solve Au=f-B^*p
1403                    #       A(u-v)=f-B^*p-Av
1404                    #       u=v+(u-v)
1405            u=v+self.solve_A(v,p)
1406    
1407              if solver=='MINRES':      
1408                    if self.verbose: print "enter MINRES method (iter_max=%s)"%max_iter
1409                    p=MINRES(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)
1410                    # solve Au=f-B^*p
1411                    #       A(u-v)=f-B^*p-Av
1412                    #       u=v+(u-v)
1413            u=v+self.solve_A(v,p)
1414                  
1415                  if solver=='PCG':
1416                    if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter
1417                    p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p)
1418                u=r[0]  
1419    
1420                  print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)
1421    
1422              return u,p
1423    
1424          def __Msolve(self,r):
1425              return self.solve_prec(r[1])
1426    
1427          def __Msolve2(self,r):
1428              return self.solve_prec(r)
1429    
1430    
1431          def __Aprod(self,p):
1432              # return BA^-1B*p
1433              #solve Av =-B^*p as Av =f-Az-B^*p
1434              v=self.solve_A(self.__z,-p)
1435              return ArithmeticTuple(v, self.B(v))
1436    
1437          def __Aprod2(self,p):
1438              # return BA^-1B*p
1439              #solve Av =-B^*p as Av =f-Az-B^*p
1440          v=self.solve_A(self.__z,-p)
1441              return self.B(v)
1442    
1443          def __Aprod_Newton(self,p):
1444              # return BA^-1B*p
1445              #solve Av =-B^*p as Av =f-Az-B^*p
1446          v=self.solve_A(self.__z,-p)
1447              return self.B(v-self.__z)
1448    
1449    class SaddlePointProblem(object):
1450       """
1451       This implements a solver for a saddlepoint problem
1452    
1453       M{f(u,p)=0}
1454       M{g(u)=0}
1455    
1456       for u and p. The problem is solved with an inexact Uszawa scheme for p:
1457    
1458       M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})}
1459       M{Q_g (p^{k+1}-p^{k}) =   g(u^{k+1})}
1460    
1461       where Q_f is an approximation of the Jacobiean A_f of f with respect to u  and Q_f is an approximation of
1462       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'
1463       Q_g can be difficult, non-linear conjugate gradient method is applied to solve for p, so Q_g plays
1464       in fact the role of a preconditioner.
1465       """
1466       def __init__(self,verbose=False,*args):
1467           """
1468           initializes the problem
1469    
1470           @param verbose: switches on the printing out some information
1471           @type verbose: C{bool}
1472           @note: this method may be overwritten by a particular saddle point problem
1473           """
1474           if not isinstance(verbose,bool):
1475                raise TypeError("verbose needs to be of type bool.")
1476           self.__verbose=verbose
1477           self.relaxation=1.
1478    
1479       def trace(self,text):
1480           """
1481           prints text if verbose has been set
1482    
1483           @param text: a text message
1484           @type text: C{str}
1485           """
1486           if self.__verbose: print "#s: #s"%(str(self),text)
1487    
1488       def solve_f(self,u,p,tol=1.e-8):
1489           """
1490           solves
1491    
1492           A_f du = f(u,p)
1493    
1494           with tolerance C{tol} and return du. A_f is Jacobiean of f with respect to u.
1495    
1496           @param u: current approximation of u
1497           @type u: L{escript.Data}
1498           @param p: current approximation of p
1499           @type p: L{escript.Data}
1500           @param tol: tolerance expected for du
1501           @type tol: C{float}
1502           @return: increment du
1503           @rtype: L{escript.Data}
1504           @note: this method has to be overwritten by a particular saddle point problem
1505           """
1506           pass
1507    
1508       def solve_g(self,u,tol=1.e-8):
1509           """
1510           solves
1511    
1512           Q_g dp = g(u)
1513    
1514           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.
1515    
1516           @param u: current approximation of u
1517           @type u: L{escript.Data}
1518           @param tol: tolerance expected for dp
1519           @type tol: C{float}
1520           @return: increment dp
1521           @rtype: L{escript.Data}
1522           @note: this method has to be overwritten by a particular saddle point problem
1523           """
1524           pass
1525    
1526       def inner(self,p0,p1):
1527           """
1528           inner product of p0 and p1 approximating p. Typically this returns integrate(p0*p1)
1529           @return: inner product of p0 and p1
1530           @rtype: C{float}
1531           """
1532           pass
1533    
1534       subiter_max=3
1535       def solve(self,u0,p0,tolerance=1.e-6,tolerance_u=None,iter_max=100,accepted_reduction=0.995,relaxation=None):
1536            """
1537            runs the solver
1538    
1539            @param u0: initial guess for C{u}
1540            @type u0: L{esys.escript.Data}
1541            @param p0: initial guess for C{p}
1542            @type p0: L{esys.escript.Data}
1543            @param tolerance: tolerance for relative error in C{u} and C{p}
1544            @type tolerance: positive C{float}
1545            @param tolerance_u: tolerance for relative error in C{u} if different from C{tolerance}
1546            @type tolerance_u: positive C{float}
1547            @param iter_max: maximum number of iteration steps.
1548            @type iter_max: C{int}
1549            @param accepted_reduction: if the norm  g cannot be reduced by C{accepted_reduction} backtracking to adapt the
1550                                       relaxation factor. If C{accepted_reduction=None} no backtracking is used.
1551            @type accepted_reduction: positive C{float} or C{None}
1552            @param relaxation: initial relaxation factor. If C{relaxation==None}, the last relaxation factor is used.
1553            @type relaxation: C{float} or C{None}
1554            """
1555            tol=1.e-2
1556            if tolerance_u==None: tolerance_u=tolerance
1557            if not relaxation==None: self.relaxation=relaxation
1558            if accepted_reduction ==None:
1559                  angle_limit=0.
1560            elif accepted_reduction>=1.:
1561                  angle_limit=0.
1562            else:
1563                  angle_limit=util.sqrt(1-accepted_reduction**2)
1564            self.iter=0
1565            u=u0
1566            p=p0
1567            #
1568            #   initialize things:
1569            #
1570            converged=False
1571            #
1572            #  start loop:
1573            #
1574            #  initial search direction is g
1575            #
1576            while not converged :
1577                if self.iter>iter_max:
1578                    raise ArithmeticError("no convergence after %s steps."%self.iter)
1579                f_new=self.solve_f(u,p,tol)
1580                norm_f_new = util.Lsup(f_new)
1581                u_new=u-f_new
1582                g_new=self.solve_g(u_new,tol)
1583                self.iter+=1
1584                norm_g_new = util.sqrt(self.inner(g_new,g_new))
1585                if norm_f_new==0. and norm_g_new==0.: return u, p
1586                if self.iter>1 and not accepted_reduction==None:
1587                   #
1588                   #   did we manage to reduce the norm of G? I
1589                   #   if not we start a backtracking procedure
1590                   #
1591                   # print "new/old norm = ",norm_g_new, norm_g, norm_g_new/norm_g
1592                   if norm_g_new > accepted_reduction * norm_g:
1593                      sub_iter=0
1594                      s=self.relaxation
1595                      d=g
1596                      g_last=g
1597                      self.trace("    start substepping: f = %s, g = %s, relaxation = %s."%(norm_f_new, norm_g_new, s))
1598                      while sub_iter < self.subiter_max and  norm_g_new > accepted_reduction * norm_g:
1599                         dg= g_new-g_last
1600                         norm_dg=abs(util.sqrt(self.inner(dg,dg))/self.relaxation)
1601                         rad=self.inner(g_new,dg)/self.relaxation
1602                         # print "   ",sub_iter,": rad, norm_dg:",abs(rad), norm_dg*norm_g_new * angle_limit
1603                         # print "   ",sub_iter,": rad, norm_dg:",rad, norm_dg, norm_g_new, norm_g
1604                         if abs(rad) < norm_dg*norm_g_new * angle_limit:
1605                             if sub_iter>0: self.trace("    no further improvements expected from backtracking.")
1606                             break
1607                         r=self.relaxation
1608                         self.relaxation= - rad/norm_dg**2
1609                         s+=self.relaxation
1610                         #####
1611                         # a=g_new+self.relaxation*dg/r
1612                         # print "predicted new norm = ",util.sqrt(self.inner(a,a)),util.sqrt(self.inner(g_new,g_new)), self.relaxation
1613                         #####
1614                         g_last=g_new
1615                         p+=self.relaxation*d
1616                         f_new=self.solve_f(u,p,tol)
1617                         u_new=u-f_new
1618                         g_new=self.solve_g(u_new,tol)
1619                         self.iter+=1
1620                         norm_f_new = util.Lsup(f_new)
1621                         norm_g_new = util.sqrt(self.inner(g_new,g_new))
1622                         # print "   ",sub_iter," new g norm",norm_g_new
1623                         self.trace("    %s th sub-step: f = %s, g = %s, relaxation = %s."%(sub_iter, norm_f_new, norm_g_new, s))
1624                         #
1625                         #   can we expect reduction of g?
1626                         #
1627                         # u_last=u_new
1628                         sub_iter+=1
1629                      self.relaxation=s
1630                #
1631                #  check for convergence:
1632                #
1633                norm_u_new = util.Lsup(u_new)
1634                p_new=p+self.relaxation*g_new
1635                norm_p_new = util.sqrt(self.inner(p_new,p_new))
1636                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))
1637    
1638                if self.iter>1:
1639                   dg2=g_new-g
1640                   df2=f_new-f
1641                   norm_dg2=util.sqrt(self.inner(dg2,dg2))
1642                   norm_df2=util.Lsup(df2)
1643                   # print norm_g_new, norm_g, norm_dg, norm_p, tolerance
1644                   tol_eq_g=tolerance*norm_dg2/(norm_g*abs(self.relaxation))*norm_p_new
1645                   tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new
1646                   if norm_g_new <= tol_eq_g and norm_f_new <= tol_eq_f:
1647                       converged=True
1648                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
1649            self.trace("convergence after %s steps."%self.iter)
1650            return u,p
1651    #   def solve(self,u0,p0,tolerance=1.e-6,iter_max=10,self.relaxation=1.):
1652    #      tol=1.e-2
1653    #      iter=0
1654    #      converged=False
1655    #      u=u0*1.
1656    #      p=p0*1.
1657    #      while not converged and iter<iter_max:
1658    #          du=self.solve_f(u,p,tol)
1659    #          u-=du
1660    #          norm_du=util.Lsup(du)
1661    #          norm_u=util.Lsup(u)
1662    #        
1663    #          dp=self.relaxation*self.solve_g(u,tol)
1664    #          p+=dp
1665    #          norm_dp=util.sqrt(self.inner(dp,dp))
1666    #          norm_p=util.sqrt(self.inner(p,p))
1667    #          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)
1668    #          iter+=1
1669    #
1670    #          converged = (norm_du <= tolerance*norm_u) and  (norm_dp <= tolerance*norm_p)
1671    #      if converged:
1672    #          print "convergence after %s steps."%iter
1673    #      else:
1674    #          raise ArithmeticError("no convergence after %s steps."%iter)
1675    #
1676    #      return u,p
1677              
1678    def MaskFromBoundaryTag(function_space,*tags):
1679       """
1680       create a mask on the given function space which one for samples
1681       that touch the boundary tagged by tags.
1682    
1683       usage: m=MaskFromBoundaryTag(Solution(domain),"left", "right")
1684    
1685       @param function_space: a given function space
1686       @type function_space: L{escript.FunctionSpace}
1687       @param tags: boundray tags
1688       @type tags: C{str}
1689       @return: a mask which marks samples used by C{function_space} that are touching the
1690                boundary tagged by any of the given tags.
1691       @rtype: L{escript.Data} of rank 0
1692       """
1693       pde=linearPDEs.LinearPDE(function_space.getDomain(),numEquations=1, numSolutions=1)
1694       d=escript.Scalar(0.,escript.FunctionOnBoundary(function_space.getDomain()))
1695       for t in tags: d.setTaggedValue(t,1.)
1696       pde.setValue(y=d)
1697       out=util.whereNonZero(pde.getRightHandSide())
1698       if out.getFunctionSpace() == function_space:
1699          return out
1700       else:
1701          return util.whereNonZero(util.interpolate(out,function_space))
1702    
1703    
1704    

Legend:
Removed from v.614  
changed lines
  Added in v.1519

  ViewVC Help
Powered by ViewVC 1.1.26