/[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 351 by gross, Tue Dec 13 09:12:15 2005 UTC revision 2061 by jfenwick, Wed Nov 19 03:40:21 2008 UTC
# Line 1  Line 1 
1  # $Id$  
2    ########################################################
3    #
4    # Copyright (c) 2003-2008 by University of Queensland
5    # Earth Systems Science Computational Center (ESSCC)
6    # http://www.uq.edu.au/esscc
7    #
8    # Primary Business: Queensland, Australia
9    # Licensed under the Open Software License version 3.0
10    # http://www.opensource.org/licenses/osl-3.0.php
11    #
12    ########################################################
13    
14    __copyright__="""Copyright (c) 2003-2008 by University of Queensland
15    Earth Systems Science Computational Center (ESSCC)
16    http://www.uq.edu.au/esscc
17    Primary Business: Queensland, Australia"""
18    __license__="""Licensed under the Open Software License version 3.0
19    http://www.opensource.org/licenses/osl-3.0.php"""
20    __url__="http://www.uq.edu.au/esscc/escript-finley"
21    
22  """  """
23  Provides some tools related to PDEs.  Provides some tools related to PDEs.
# Line 7  Currently includes: Line 26  Currently includes:
26      - Projector - to project a discontinuous      - Projector - to project a discontinuous
27      - Locator - to trace values in data objects at a certain location      - Locator - to trace values in data objects at a certain location
28      - TimeIntegrationManager - to handel extraplotion in time      - TimeIntegrationManager - to handel extraplotion in time
29        - SaddlePointProblem - solver for Saddle point problems using the inexact uszawa scheme
30    
31    @var __author__: name of author
32    @var __copyright__: copyrights
33    @var __license__: licence agreement
34    @var __url__: url entry point on documentation
35    @var __version__: version
36    @var __date__: date of the version
37  """  """
38    
39    __author__="Lutz Gross, l.gross@uq.edu.au"
40    
41    
42  import escript  import escript
43  import linearPDEs  import linearPDEs
44  import numarray  import numarray
45  import util  import util
46    import math
47    
48    ##### Added by Artak
49    # from Numeric import zeros,Int,Float64
50    ###################################
51    
52    
53  class TimeIntegrationManager:  class TimeIntegrationManager:
54    """    """
55    a simple mechanism to manage time dependend values.    a simple mechanism to manage time dependend values.
56    
57    typical usage is:    typical usage is::
58    
59    dt=0.1 # time increment       dt=0.1 # time increment
60    tm=TimeIntegrationManager(inital_value,p=1)       tm=TimeIntegrationManager(inital_value,p=1)
61    while t<1.       while t<1.
62        v_guess=tm.extrapolate(dt) # extrapolate to t+dt           v_guess=tm.extrapolate(dt) # extrapolate to t+dt
63        v=...           v=...
64        tm.checkin(dt,v)           tm.checkin(dt,v)
65        t+=dt           t+=dt
66    
67    @remark: currently only p=1 is supported.    @note: currently only p=1 is supported.
68    """    """
69    def __init__(self,*inital_values,**kwargs):    def __init__(self,*inital_values,**kwargs):
70       """       """
# Line 49  class TimeIntegrationManager: Line 85  class TimeIntegrationManager:
85    
86    def getTime(self):    def getTime(self):
87        return self.__t        return self.__t
88      def getValue(self):
89          out=self.__v_mem[0]
90          if len(out)==1:
91              return out[0]
92          else:
93              return out
94    
95    def checkin(self,dt,*values):    def checkin(self,dt,*values):
96        """        """
# Line 84  class TimeIntegrationManager: Line 126  class TimeIntegrationManager:
126           return out[0]           return out[0]
127        else:        else:
128           return out           return out
129    
130    
131  class Projector:  class Projector:
132    """    """
# Line 106  class Projector: Line 149  class Projector:
149      self.__pde.setValue(D = 1.)      self.__pde.setValue(D = 1.)
150      return      return
151    
   def __del__(self):  
     return  
   
152    def __call__(self, input_data):    def __call__(self, input_data):
153      """      """
154      Projects input_data onto a continuous function      Projects input_data onto a continuous function
155    
156      @param input_data: The input_data to be projected.      @param input_data: The input_data to be projected.
157      """      """
158      out=escript.Data(0.,input_data.getShape(),what=escript.ContinuousFunction(self.__pde.getDomain()))      out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())
159        self.__pde.setValue(Y = escript.Data(), Y_reduced = escript.Data())
160      if input_data.getRank()==0:      if input_data.getRank()==0:
161          self.__pde.setValue(Y = input_data)          self.__pde.setValue(Y = input_data)
162          out=self.__pde.getSolution()          out=self.__pde.getSolution()
# Line 143  class Projector: Line 184  class Projector:
184                      out[i0,i1,i2,i3]=self.__pde.getSolution()                      out[i0,i1,i2,i3]=self.__pde.getSolution()
185      return out      return out
186    
187    class NoPDE:
188         """
189         solves the following problem for u:
190    
191         M{kronecker[i,j]*D[j]*u[j]=Y[i]}
192    
193         with constraint
194    
195         M{u[j]=r[j]}  where M{q[j]>0}
196    
197         where D, Y, r and q are given functions of rank 1.
198    
199         In the case of scalars this takes the form
200    
201         M{D*u=Y}
202    
203         with constraint
204    
205         M{u=r}  where M{q>0}
206    
207         where D, Y, r and q are given scalar functions.
208    
209         The constraint is overwriting any other condition.
210    
211         @note: This class is similar to the L{linearPDEs.LinearPDE} class with A=B=C=X=0 but has the intention
212                that all input parameter are given in L{Solution} or L{ReducedSolution}. The whole
213                thing is a bit strange and I blame Robert.Woodcock@csiro.au for this.
214         """
215         def __init__(self,domain,D=None,Y=None,q=None,r=None):
216             """
217             initialize the problem
218    
219             @param domain: domain of the PDE.
220             @type domain: L{Domain}
221             @param D: coefficient of the solution.
222             @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}
223             @param Y: right hand side
224             @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}
225             @param q: location of constraints
226             @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}
227             @param r: value of solution at locations of constraints
228             @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}
229             """
230             self.__domain=domain
231             self.__D=D
232             self.__Y=Y
233             self.__q=q
234             self.__r=r
235             self.__u=None
236             self.__function_space=escript.Solution(self.__domain)
237         def setReducedOn(self):
238             """
239             sets the L{FunctionSpace} of the solution to L{ReducedSolution}
240             """
241             self.__function_space=escript.ReducedSolution(self.__domain)
242             self.__u=None
243    
244         def setReducedOff(self):
245             """
246             sets the L{FunctionSpace} of the solution to L{Solution}
247             """
248             self.__function_space=escript.Solution(self.__domain)
249             self.__u=None
250            
251         def setValue(self,D=None,Y=None,q=None,r=None):
252             """
253             assigns values to the parameters.
254    
255             @param D: coefficient of the solution.
256             @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}
257             @param Y: right hand side
258             @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}
259             @param q: location of constraints
260             @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}
261             @param r: value of solution at locations of constraints
262             @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}
263             """
264             if not D==None:
265                self.__D=D
266                self.__u=None
267             if not Y==None:
268                self.__Y=Y
269                self.__u=None
270             if not q==None:
271                self.__q=q
272                self.__u=None
273             if not r==None:
274                self.__r=r
275                self.__u=None
276    
277         def getSolution(self):
278             """
279             returns the solution
280            
281             @return: the solution of the problem
282             @rtype: L{Data} object in the L{FunctionSpace} L{Solution} or L{ReducedSolution}.
283             """
284             if self.__u==None:
285                if self.__D==None:
286                   raise ValueError,"coefficient D is undefined"
287                D=escript.Data(self.__D,self.__function_space)
288                if D.getRank()>1:
289                   raise ValueError,"coefficient D must have rank 0 or 1"
290                if self.__Y==None:
291                   self.__u=escript.Data(0.,D.getShape(),self.__function_space)
292                else:
293                   self.__u=util.quotient(self.__Y,D)
294                if not self.__q==None:
295                    q=util.wherePositive(escript.Data(self.__q,self.__function_space))
296                    self.__u*=(1.-q)
297                    if not self.__r==None: self.__u+=q*self.__r
298             return self.__u
299                
300  class Locator:  class Locator:
301       """       """
302       Locator provides access to the values of data objects at a given       Locator provides access to the values of data objects at a given
# Line 159  class Locator: Line 312  class Locator:
312         Initializes a Locator to access values in Data objects on the Doamin         Initializes a Locator to access values in Data objects on the Doamin
313         or FunctionSpace where for the sample point which         or FunctionSpace where for the sample point which
314         closest to the given point x.         closest to the given point x.
315    
316           @param where: function space
317           @type where: L{escript.FunctionSpace}
318           @param x: coefficient of the solution.
319           @type x: L{numarray.NumArray} or C{list} of L{numarray.NumArray}
320         """         """
321         if isinstance(where,escript.FunctionSpace):         if isinstance(where,escript.FunctionSpace):
322            self.__function_space=where            self.__function_space=where
323         else:         else:
324            self.__function_space=escript.ContinuousFunction(where)            self.__function_space=escript.ContinuousFunction(where)
325         self.__id=util.length(x[:self.__function_space.getDim()]-self.__function_space.getX()).mindp()         if isinstance(x, list):
326               self.__id=[]
327               for p in x:
328                  self.__id.append(util.length(self.__function_space.getX()-p[:self.__function_space.getDim()]).minGlobalDataPoint())
329           else:
330               self.__id=util.length(self.__function_space.getX()-x[:self.__function_space.getDim()]).minGlobalDataPoint()
331    
332       def __str__(self):       def __str__(self):
333         """         """
334         Returns the coordinates of the Locator as a string.         Returns the coordinates of the Locator as a string.
335         """         """
336         return "<Locator %s>"%str(self.getX())         x=self.getX()
337           if instance(x,list):
338              out="["
339              first=True
340              for xx in x:
341                if not first:
342                    out+=","
343                else:
344                    first=False
345                out+=str(xx)
346              out+="]>"
347           else:
348              out=str(x)
349           return out
350    
351         def getX(self):
352            """
353        Returns the exact coordinates of the Locator.
354        """
355            return self(self.getFunctionSpace().getX())
356    
357       def getFunctionSpace(self):       def getFunctionSpace(self):
358          """          """
# Line 178  class Locator: Line 360  class Locator:
360      """      """
361          return self.__function_space          return self.__function_space
362    
363       def getId(self):       def getId(self,item=None):
364          """          """
365      Returns the identifier of the location.      Returns the identifier of the location.
366      """      """
367          return self.__id          if item == None:
368               return self.__id
369            else:
370               if isinstance(self.__id,list):
371                  return self.__id[item]
372               else:
373                  return self.__id
374    
      def getX(self):  
         """  
     Returns the exact coordinates of the Locator.  
     """  
         return self(self.getFunctionSpace().getX())  
375    
376       def __call__(self,data):       def __call__(self,data):
377          """          """
# Line 204  class Locator: Line 387  class Locator:
387      """      """
388          if isinstance(data,escript.Data):          if isinstance(data,escript.Data):
389             if data.getFunctionSpace()==self.getFunctionSpace():             if data.getFunctionSpace()==self.getFunctionSpace():
390               out=data.convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])               dat=data
391             else:             else:
392               out=data.interpolate(self.getFunctionSpace()).convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])               dat=data.interpolate(self.getFunctionSpace())
393             if data.getRank()==0:             id=self.getId()
394                return out[0]             r=data.getRank()
395               if isinstance(id,list):
396                   out=[]
397                   for i in id:
398                      o=data.getValueOfGlobalDataPoint(*i)
399                      if data.getRank()==0:
400                         out.append(o[0])
401                      else:
402                         out.append(o)
403                   return out
404             else:             else:
405                return out               out=data.getValueOfGlobalDataPoint(*id)
406                 if data.getRank()==0:
407                    return out[0]
408                 else:
409                    return out
410          else:          else:
411             return data             return data
412    
413  # vim: expandtab shiftwidth=4:  class SolverSchemeException(Exception):
414       """
415       exceptions thrown by solvers
416       """
417       pass
418    
419    class IndefinitePreconditioner(SolverSchemeException):
420       """
421       the preconditioner is not positive definite.
422       """
423       pass
424    class MaxIterReached(SolverSchemeException):
425       """
426       maxium number of iteration steps is reached.
427       """
428       pass
429    class IterationBreakDown(SolverSchemeException):
430       """
431       iteration scheme econouters an incurable breakdown.
432       """
433       pass
434    class NegativeNorm(SolverSchemeException):
435       """
436       a norm calculation returns a negative norm.
437       """
438       pass
439    
440    class IterationHistory(object):
441       """
442       The IterationHistory class is used to define a stopping criterium. It keeps track of the
443       residual norms. The stoppingcriterium indicates termination if the residual norm has been reduced by
444       a given tolerance.
445       """
446       def __init__(self,tolerance=math.sqrt(util.EPSILON),verbose=False):
447          """
448          Initialization
449    
450          @param tolerance: tolerance
451          @type tolerance: positive C{float}
452          @param verbose: switches on the printing out some information
453          @type verbose: C{bool}
454          """
455          if not tolerance>0.:
456              raise ValueError,"tolerance needs to be positive."
457          self.tolerance=tolerance
458          self.verbose=verbose
459          self.history=[]
460       def stoppingcriterium(self,norm_r,r,x):
461           """
462           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.
463    
464          
465           @param norm_r: current residual norm
466           @type norm_r: non-negative C{float}
467           @param r: current residual (not used)
468           @param x: current solution approximation (not used)
469           @return: C{True} is the stopping criterium is fullfilled. Otherwise C{False} is returned.
470           @rtype: C{bool}
471    
472           """
473           self.history.append(norm_r)
474           if self.verbose: print "iter: %s:  inner(rhat,r) = %e"#(len(self.history)-1, self.history[-1])
475           return self.history[-1]<=self.tolerance * self.history[0]
476    
477       def stoppingcriterium2(self,norm_r,norm_b,solver="GMRES",TOL=None):
478           """
479           returns True if the C{norm_r} is C{tolerance}*C{norm_b}
480    
481          
482           @param norm_r: current residual norm
483           @type norm_r: non-negative C{float}
484           @param norm_b: norm of right hand side
485           @type norm_b: non-negative C{float}
486           @return: C{True} is the stopping criterium is fullfilled. Otherwise C{False} is returned.
487           @rtype: C{bool}
488    
489           """
490           if TOL==None:
491              TOL=self.tolerance
492           self.history.append(norm_r)
493           if self.verbose: print "iter: %s:  norm(r) = %e"#(len(self.history)-1, self.history[-1])
494           return self.history[-1]<=TOL * norm_b
495    
496    def PCG(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
497       """
498       Solver for
499    
500       M{Ax=b}
501    
502       with a symmetric and positive definite operator A (more details required!).
503       It uses the conjugate gradient method with preconditioner M providing an approximation of A.
504    
505       The iteration is terminated if the C{stoppingcriterium} function return C{True}.
506    
507       For details on the preconditioned conjugate gradient method see the book:
508    
509       Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
510       T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
511       C. Romine, and H. van der Vorst.
512    
513       @param b: the right hand side of the liner system. C{b} is altered.
514       @type b: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
515       @param Aprod: returns the value Ax
516       @type Aprod: function C{Aprod(x)} where C{x} is of the same object like argument C{x}.
517           The returned object needs to be of the same type like argument C{b}.
518       @param Msolve: solves Mx=r
519       @type Msolve: function C{Msolve(r)} where C{r} is of the same type like argument C{b}.
520          The returned object needs to be of the same type like argument C{x}.
521       @param bilinearform: inner product C{<x,r>}
522       @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same type like argument C{x} and C{r} is.
523           The returned value is a C{float}.
524       @param stoppingcriterium: function which returns True if a stopping criterium is meet.
525           C{stoppingcriterium} has the arguments C{norm_r}, C{r} and C{x} giving the current
526           norm of the residual (=C{sqrt(bilinearform(Msolve(r),r)}), the current residual and
527           the current solution approximation. C{stoppingcriterium} is called in each iteration step.
528       @type stoppingcriterium: function that returns C{True} or C{False}
529       @param x: an initial guess for the solution. If no C{x} is given 0*b is used.
530       @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
531       @param iter_max: maximum number of iteration steps.
532       @type iter_max: C{int}
533       @return: the solution approximation and the corresponding residual
534       @rtype: C{tuple}
535       @warning: C{b} and C{x} are altered.
536       """
537       iter=0
538       if x==None:
539          x=0*b
540       else:
541          b += (-1)*Aprod(x)
542       r=b
543       rhat=Msolve(r)
544       d = rhat
545       rhat_dot_r = bilinearform(rhat, r)
546       if rhat_dot_r<0: raise NegativeNorm,"negative norm."
547    
548       while not stoppingcriterium(math.sqrt(rhat_dot_r),r,x):
549           iter+=1
550           if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
551    
552           q=Aprod(d)
553           alpha = rhat_dot_r / bilinearform(d, q)
554           x += alpha * d
555           r += (-alpha) * q
556    
557           rhat=Msolve(r)
558           rhat_dot_r_new = bilinearform(rhat, r)
559           beta = rhat_dot_r_new / rhat_dot_r
560           rhat+=beta * d
561           d=rhat
562    
563           rhat_dot_r = rhat_dot_r_new
564           if rhat_dot_r<0: raise NegativeNorm,"negative norm."
565    
566       return x,r
567    
568    class Defect(object):
569        """
570        defines a non-linear defect F(x) of a variable x
571        """
572        def __init__(self):
573            """
574            initialize defect
575            """
576            self.setDerivativeIncrementLength()
577    
578        def bilinearform(self, x0, x1):
579            """
580            returns the inner product of x0 and x1
581            @param x0: a value for x
582            @param x1: a value for x
583            @return: the inner product of x0 and x1
584            @rtype: C{float}
585            """
586            return 0
587          
588        def norm(self,x):
589            """
590            the norm of argument C{x}
591    
592            @param x: a value for x
593            @return: norm of argument x
594            @rtype: C{float}
595            @note: by default C{sqrt(self.bilinearform(x,x)} is retrurned.
596            """
597            s=self.bilinearform(x,x)
598            if s<0: raise NegativeNorm,"negative norm."
599            return math.sqrt(s)
600    
601    
602        def eval(self,x):
603            """
604            returns the value F of a given x
605    
606            @param x: value for which the defect C{F} is evalulated.
607            @return: value of the defect at C{x}
608            """
609            return 0
610    
611        def __call__(self,x):
612            return self.eval(x)
613    
614        def setDerivativeIncrementLength(self,inc=math.sqrt(util.EPSILON)):
615            """
616            sets the relative length of the increment used to approximate the derivative of the defect
617            the increment is inc*norm(x)/norm(v)*v in the direction of v with x as a staring point.
618    
619            @param inc: relative increment length
620            @type inc: positive C{float}
621            """
622            if inc<=0: raise ValueError,"positive increment required."
623            self.__inc=inc
624    
625        def getDerivativeIncrementLength(self):
626            """
627            returns the relative increment length used to approximate the derivative of the defect
628            @return: value of the defect at C{x}
629            @rtype: positive C{float}
630            """
631            return self.__inc
632    
633        def derivative(self, F0, x0, v, v_is_normalised=True):
634            """
635            returns the directional derivative at x0 in the direction of v
636    
637            @param F0: value of this defect at x0
638            @param x0: value at which derivative is calculated.
639            @param v: direction
640            @param v_is_normalised: is true to indicate that C{v} is nomalized (self.norm(v)=0)
641            @return: derivative of this defect at x0 in the direction of C{v}
642            @note: by default numerical evaluation (self.eval(x0+eps*v)-F0)/eps is used but this method
643            maybe oepsnew verwritten to use exact evalution.
644            """
645            normx=self.norm(x0)
646            if normx>0:
647                 epsnew = self.getDerivativeIncrementLength() * normx
648            else:
649                 epsnew = self.getDerivativeIncrementLength()
650            if not v_is_normalised:
651                normv=self.norm(v)
652                if normv<=0:
653                   return F0*0
654                else:
655                   epsnew /= normv
656            F1=self.eval(x0 + epsnew * v)
657            return (F1-F0)/epsnew
658    
659    ######################################    
660    def NewtonGMRES(defect, x, iter_max=100, sub_iter_max=20, atol=0,rtol=1.e-4, sub_tol_max=0.5, gamma=0.9, verbose=False):
661       """
662       solves a non-linear problem M{F(x)=0} for unknown M{x} using the stopping criterion:
663    
664       M{norm(F(x) <= atol + rtol * norm(F(x0)}
665      
666       where M{x0} is the initial guess.
667    
668       @param defect: object defining the the function M{F}, C{defect.norm} defines the M{norm} used in the stopping criterion.
669       @type defect: L{Defect}
670       @param x: initial guess for the solution, C{x} is altered.
671       @type x: any object type allowing basic operations such as  L{numarray.NumArray}, L{Data}
672       @param iter_max: maximum number of iteration steps
673       @type iter_max: positive C{int}
674       @param sub_iter_max:
675       @type sub_iter_max:
676       @param atol: absolute tolerance for the solution
677       @type atol: positive C{float}
678       @param rtol: relative tolerance for the solution
679       @type rtol: positive C{float}
680       @param gamma: tolerance safety factor for inner iteration
681       @type gamma: positive C{float}, less than 1
682       @param sub_tol_max: upper bound for inner tolerance.
683       @type sub_tol_max: positive C{float}, less than 1
684       @return: an approximation of the solution with the desired accuracy
685       @rtype: same type as the initial guess.
686       """
687       lmaxit=iter_max
688       if atol<0: raise ValueError,"atol needs to be non-negative."
689       if rtol<0: raise ValueError,"rtol needs to be non-negative."
690       if rtol+atol<=0: raise ValueError,"rtol or atol needs to be non-negative."
691       if gamma<=0 or gamma>=1: raise ValueError,"tolerance safety factor for inner iteration (gamma =%s) needs to be positive and less than 1."%gamma
692       if sub_tol_max<=0 or sub_tol_max>=1: raise ValueError,"upper bound for inner tolerance for inner iteration (sub_tol_max =%s) needs to be positive and less than 1."%sub_tol_max
693    
694       F=defect(x)
695       fnrm=defect.norm(F)
696       stop_tol=atol + rtol*fnrm
697       sub_tol=sub_tol_max
698       if verbose: print "NewtonGMRES: initial residual = %e."%fnrm
699       if verbose: print "             tolerance = %e."%sub_tol
700       iter=1
701       #
702       # main iteration loop
703       #
704       while not fnrm<=stop_tol:
705                if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
706                #
707            #   adjust sub_tol_
708            #
709                if iter > 1:
710               rat=fnrm/fnrmo
711                   sub_tol_old=sub_tol
712               sub_tol=gamma*rat**2
713               if gamma*sub_tol_old**2 > .1: sub_tol=max(sub_tol,gamma*sub_tol_old**2)
714               sub_tol=max(min(sub_tol,sub_tol_max), .5*stop_tol/fnrm)
715            #
716            # calculate newton increment xc
717                #     if iter_max in __FDGMRES is reached MaxIterReached is thrown
718                #     if iter_restart -1 is returned as sub_iter
719                #     if  atol is reached sub_iter returns the numer of steps performed to get there
720                #
721                #  
722                if verbose: print "             subiteration (GMRES) is called with relative tolerance %e."%sub_tol
723                try:
724                   xc, sub_iter=__FDGMRES(F, defect, x, sub_tol*fnrm, iter_max=iter_max-iter, iter_restart=sub_iter_max)
725                except MaxIterReached:
726                   raise MaxIterReached,"maximum number of %s steps reached."%iter_max
727                if sub_iter<0:
728                   iter+=sub_iter_max
729                else:
730                   iter+=sub_iter
731                # ====
732            x+=xc
733                F=defect(x)
734            iter+=1
735                fnrmo, fnrm=fnrm, defect.norm(F)
736                if verbose: print "             step %s: residual %e."%(iter,fnrm)
737       if verbose: print "NewtonGMRES: completed after %s steps."%iter
738       return x
739    
740    def __givapp(c,s,vin):
741        """
742        apply a sequence of Givens rotations (c,s) to the recuirsively to the vector vin
743        @warning: C{vin} is altered.
744        """
745        vrot=vin
746        if isinstance(c,float):
747            vrot=[c*vrot[0]-s*vrot[1],s*vrot[0]+c*vrot[1]]
748        else:
749            for i in range(len(c)):
750                w1=c[i]*vrot[i]-s[i]*vrot[i+1]
751            w2=s[i]*vrot[i]+c[i]*vrot[i+1]
752                vrot[i:i+2]=w1,w2
753        return vrot
754    
755    def __FDGMRES(F0, defect, x0, atol, iter_max=100, iter_restart=20):
756       h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)
757       c=numarray.zeros(iter_restart,numarray.Float64)
758       s=numarray.zeros(iter_restart,numarray.Float64)
759       g=numarray.zeros(iter_restart,numarray.Float64)
760       v=[]
761    
762       rho=defect.norm(F0)
763       if rho<=0.: return x0*0
764      
765       v.append(-F0/rho)
766       g[0]=rho
767       iter=0
768       while rho > atol and iter<iter_restart-1:
769    
770        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
771    
772            p=defect.derivative(F0,x0,v[iter], v_is_normalised=True)
773        v.append(p)
774    
775        v_norm1=defect.norm(v[iter+1])
776    
777            # Modified Gram-Schmidt
778        for j in range(iter+1):
779             h[j][iter]=defect.bilinearform(v[j],v[iter+1])  
780             v[iter+1]-=h[j][iter]*v[j]
781          
782        h[iter+1][iter]=defect.norm(v[iter+1])
783        v_norm2=h[iter+1][iter]
784    
785            # Reorthogonalize if needed
786        if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
787            for j in range(iter+1):  
788               hr=defect.bilinearform(v[j],v[iter+1])
789                   h[j][iter]=h[j][iter]+hr
790                   v[iter+1] -= hr*v[j]
791    
792            v_norm2=defect.norm(v[iter+1])
793            h[iter+1][iter]=v_norm2
794            #   watch out for happy breakdown
795            if not v_norm2 == 0:
796                    v[iter+1]=v[iter+1]/h[iter+1][iter]
797    
798            #   Form and store the information for the new Givens rotation
799        if iter > 0 :
800            hhat=numarray.zeros(iter+1,numarray.Float64)
801            for i in range(iter+1) : hhat[i]=h[i][iter]
802            hhat=__givapp(c[0:iter],s[0:iter],hhat);
803                for i in range(iter+1) : h[i][iter]=hhat[i]
804    
805        mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])
806    
807        if mu!=0 :
808            c[iter]=h[iter][iter]/mu
809            s[iter]=-h[iter+1][iter]/mu
810            h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter]
811            h[iter+1][iter]=0.0
812            g[iter:iter+2]=__givapp(c[iter],s[iter],g[iter:iter+2])
813    
814            # Update the residual norm
815            rho=abs(g[iter+1])
816        iter+=1
817    
818       # At this point either iter > iter_max or rho < tol.
819       # It's time to compute x and leave.        
820       if iter > 0 :
821         y=numarray.zeros(iter,numarray.Float64)    
822         y[iter-1] = g[iter-1] / h[iter-1][iter-1]
823         if iter > 1 :  
824            i=iter-2  
825            while i>=0 :
826              y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]
827              i=i-1
828         xhat=v[iter-1]*y[iter-1]
829         for i in range(iter-1):
830        xhat += v[i]*y[i]
831       else :
832          xhat=v[0] * 0
833    
834       if iter<iter_restart-1:
835          stopped=iter
836       else:
837          stopped=-1
838    
839       return xhat,stopped
840    
841    ##############################################
842    def GMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):
843    ################################################
844       m=iter_restart
845       iter=0
846       xc=x
847       while True:
848          if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached"%iter_max
849          xc,stopped=__GMRESm(b*1, Aprod, Msolve, bilinearform, stoppingcriterium, x=xc*1, iter_max=iter_max-iter, iter_restart=m)
850          if stopped: break
851          iter+=iter_restart    
852       return xc
853    
854    def __GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):
855       iter=0
856       r=Msolve(b)
857       r_dot_r = bilinearform(r, r)
858       if r_dot_r<0: raise NegativeNorm,"negative norm."
859       norm_b=math.sqrt(r_dot_r)
860    
861       if x==None:
862          x=0*b
863       else:
864          r=Msolve(b-Aprod(x))
865          r_dot_r = bilinearform(r, r)
866          if r_dot_r<0: raise NegativeNorm,"negative norm."
867      
868       h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)
869       c=numarray.zeros(iter_restart,numarray.Float64)
870       s=numarray.zeros(iter_restart,numarray.Float64)
871       g=numarray.zeros(iter_restart,numarray.Float64)
872       v=[]
873    
874       rho=math.sqrt(r_dot_r)
875      
876       v.append(r/rho)
877       g[0]=rho
878    
879       while not (stoppingcriterium(rho,norm_b) or iter==iter_restart-1):
880    
881        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
882    
883        p=Msolve(Aprod(v[iter]))
884    
885        v.append(p)
886    
887        v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
888    
889    # Modified Gram-Schmidt
890        for j in range(iter+1):
891          h[j][iter]=bilinearform(v[j],v[iter+1])  
892          v[iter+1]-=h[j][iter]*v[j]
893          
894        h[iter+1][iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
895        v_norm2=h[iter+1][iter]
896    
897    # Reorthogonalize if needed
898        if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
899         for j in range(iter+1):  
900            hr=bilinearform(v[j],v[iter+1])
901                h[j][iter]=h[j][iter]+hr
902                v[iter+1] -= hr*v[j]
903    
904         v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
905         h[iter+1][iter]=v_norm2
906    
907    #   watch out for happy breakdown
908            if not v_norm2 == 0:
909             v[iter+1]=v[iter+1]/h[iter+1][iter]
910    
911    #   Form and store the information for the new Givens rotation
912        if iter > 0 :
913            hhat=numarray.zeros(iter+1,numarray.Float64)
914            for i in range(iter+1) : hhat[i]=h[i][iter]
915            hhat=__givapp(c[0:iter],s[0:iter],hhat);
916                for i in range(iter+1) : h[i][iter]=hhat[i]
917    
918        mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])
919    
920        if mu!=0 :
921            c[iter]=h[iter][iter]/mu
922            s[iter]=-h[iter+1][iter]/mu
923            h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter]
924            h[iter+1][iter]=0.0
925            g[iter:iter+2]=__givapp(c[iter],s[iter],g[iter:iter+2])
926    
927    # Update the residual norm
928                  
929            rho=abs(g[iter+1])
930        iter+=1
931    
932    # At this point either iter > iter_max or rho < tol.
933    # It's time to compute x and leave.        
934    
935       if iter > 0 :
936         y=numarray.zeros(iter,numarray.Float64)    
937         y[iter-1] = g[iter-1] / h[iter-1][iter-1]
938         if iter > 1 :  
939            i=iter-2  
940            while i>=0 :
941              y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]
942              i=i-1
943         xhat=v[iter-1]*y[iter-1]
944         for i in range(iter-1):
945        xhat += v[i]*y[i]
946       else : xhat=v[0]
947    
948       x += xhat
949       if iter<iter_restart-1:
950          stopped=True
951       else:
952          stopped=False
953    
954       return x,stopped
955    
956    #################################################
957    def MINRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
958    #################################################
959        #
960        #  minres solves the system of linear equations Ax = b
961        #  where A is a symmetric matrix (possibly indefinite or singular)
962        #  and b is a given vector.
963        #  
964        #  "A" may be a dense or sparse matrix (preferably sparse!)
965        #  or the name of a function such that
966        #               y = A(x)
967        #  returns the product y = Ax for any given vector x.
968        #
969        #  "M" defines a positive-definite preconditioner M = C C'.
970        #  "M" may be a dense or sparse matrix (preferably sparse!)
971        #  or the name of a function such that
972        #  solves the system My = x for any given vector x.
973        #
974        #
975        
976        #------------------------------------------------------------------
977        # Set up y and v for the first Lanczos vector v1.
978        # y  =  beta1 P' v1,  where  P = C**(-1).
979        # v is really P' v1.
980        #------------------------------------------------------------------
981        if x==None:
982          x=0*b
983        else:
984          b += (-1)*Aprod(x)
985    
986        r1    = b
987        y = Msolve(b)
988        beta1 = bilinearform(y,b)
989    
990        if beta1< 0: raise NegativeNorm,"negative norm."
991    
992        #  If b = 0 exactly, stop with x = 0.
993        if beta1==0: return x*0.
994    
995        if beta1> 0:
996          beta1  = math.sqrt(beta1)      
997    
998        #------------------------------------------------------------------
999        # Initialize quantities.
1000        # ------------------------------------------------------------------
1001        iter   = 0
1002        Anorm = 0
1003        ynorm = 0
1004        oldb   = 0
1005        beta   = beta1
1006        dbar   = 0
1007        epsln  = 0
1008        phibar = beta1
1009        rhs1   = beta1
1010        rhs2   = 0
1011        rnorm  = phibar
1012        tnorm2 = 0
1013        ynorm2 = 0
1014        cs     = -1
1015        sn     = 0
1016        w      = b*0.
1017        w2     = b*0.
1018        r2     = r1
1019        eps    = 0.0001
1020    
1021        #---------------------------------------------------------------------
1022        # Main iteration loop.
1023        # --------------------------------------------------------------------
1024        while not stoppingcriterium(rnorm,Anorm*ynorm,'MINRES'):    #  checks ||r|| < (||A|| ||x||) * TOL
1025    
1026        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
1027            iter    = iter  +  1
1028    
1029            #-----------------------------------------------------------------
1030            # Obtain quantities for the next Lanczos vector vk+1, k = 1, 2,...
1031            # The general iteration is similar to the case k = 1 with v0 = 0:
1032            #
1033            #   p1      = Operator * v1  -  beta1 * v0,
1034            #   alpha1  = v1'p1,
1035            #   q2      = p2  -  alpha1 * v1,
1036            #   beta2^2 = q2'q2,
1037            #   v2      = (1/beta2) q2.
1038            #
1039            # Again, y = betak P vk,  where  P = C**(-1).
1040            #-----------------------------------------------------------------
1041            s = 1/beta                 # Normalize previous vector (in y).
1042            v = s*y                    # v = vk if P = I
1043        
1044            y      = Aprod(v)
1045        
1046            if iter >= 2:
1047              y = y - (beta/oldb)*r1
1048    
1049            alfa   = bilinearform(v,y)              # alphak
1050            y      += (- alfa/beta)*r2
1051            r1     = r2
1052            r2     = y
1053            y = Msolve(r2)
1054            oldb   = beta                           # oldb = betak
1055            beta   = bilinearform(y,r2)             # beta = betak+1^2
1056            if beta < 0: raise NegativeNorm,"negative norm."
1057    
1058            beta   = math.sqrt( beta )
1059            tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta
1060            
1061            if iter==1:                 # Initialize a few things.
1062              gmax   = abs( alfa )      # alpha1
1063              gmin   = gmax             # alpha1
1064    
1065            # Apply previous rotation Qk-1 to get
1066            #   [deltak epslnk+1] = [cs  sn][dbark    0   ]
1067            #   [gbar k dbar k+1]   [sn -cs][alfak betak+1].
1068        
1069            oldeps = epsln
1070            delta  = cs * dbar  +  sn * alfa  # delta1 = 0         deltak
1071            gbar   = sn * dbar  -  cs * alfa  # gbar 1 = alfa1     gbar k
1072            epsln  =               sn * beta  # epsln2 = 0         epslnk+1
1073            dbar   =            -  cs * beta  # dbar 2 = beta2     dbar k+1
1074    
1075            # Compute the next plane rotation Qk
1076    
1077            gamma  = math.sqrt(gbar*gbar+beta*beta)  # gammak
1078            gamma  = max(gamma,eps)
1079            cs     = gbar / gamma             # ck
1080            sn     = beta / gamma             # sk
1081            phi    = cs * phibar              # phik
1082            phibar = sn * phibar              # phibark+1
1083    
1084            # Update  x.
1085    
1086            denom = 1/gamma
1087            w1    = w2
1088            w2    = w
1089            w     = (v - oldeps*w1 - delta*w2) * denom
1090            x     +=  phi*w
1091    
1092            # Go round again.
1093    
1094            gmax   = max(gmax,gamma)
1095            gmin   = min(gmin,gamma)
1096            z      = rhs1 / gamma
1097            ynorm2 = z*z  +  ynorm2
1098            rhs1   = rhs2 -  delta*z
1099            rhs2   =      -  epsln*z
1100    
1101            # Estimate various norms and test for convergence.
1102    
1103            Anorm  = math.sqrt( tnorm2 )
1104            ynorm  = math.sqrt( ynorm2 )
1105    
1106            rnorm  = phibar
1107    
1108        return x
1109    
1110    def TFQMR(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
1111    
1112    # TFQMR solver for linear systems
1113    #
1114    #
1115    # initialization
1116    #
1117      errtol = math.sqrt(bilinearform(b,b))
1118      norm_b=errtol
1119      kmax  = iter_max
1120      error = []
1121    
1122      if math.sqrt(bilinearform(x,x)) != 0.0:
1123        r = b - Aprod(x)
1124      else:
1125        r = b
1126    
1127      r=Msolve(r)
1128    
1129      u1=0
1130      u2=0
1131      y1=0
1132      y2=0
1133    
1134      w = r
1135      y1 = r
1136      iter = 0
1137      d = 0
1138      
1139      v = Msolve(Aprod(y1))
1140      u1 = v
1141      
1142      theta = 0.0;
1143      eta = 0.0;
1144      tau = math.sqrt(bilinearform(r,r))
1145      error = [ error, tau ]
1146      rho = tau * tau
1147      m=1
1148    #
1149    #  TFQMR iteration
1150    #
1151    #  while ( iter < kmax-1 ):
1152      
1153      while not stoppingcriterium(tau*math.sqrt ( m + 1 ),norm_b,'TFQMR'):
1154        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
1155    
1156        sigma = bilinearform(r,v)
1157    
1158        if ( sigma == 0.0 ):
1159          raise 'TFQMR breakdown, sigma=0'
1160        
1161    
1162        alpha = rho / sigma
1163    
1164        for j in range(2):
1165    #
1166    #   Compute y2 and u2 only if you have to
1167    #
1168          if ( j == 1 ):
1169            y2 = y1 - alpha * v
1170            u2 = Msolve(Aprod(y2))
1171    
1172          m = 2 * (iter+1) - 2 + (j+1)
1173          if j==0:
1174             w = w - alpha * u1
1175             d = y1 + ( theta * theta * eta / alpha ) * d
1176          if j==1:
1177             w = w - alpha * u2
1178             d = y2 + ( theta * theta * eta / alpha ) * d
1179    
1180          theta = math.sqrt(bilinearform(w,w))/ tau
1181          c = 1.0 / math.sqrt ( 1.0 + theta * theta )
1182          tau = tau * theta * c
1183          eta = c * c * alpha
1184          x = x + eta * d
1185    #
1186    #  Try to terminate the iteration at each pass through the loop
1187    #
1188         # if ( tau * math.sqrt ( m + 1 ) <= errtol ):
1189         #   error = [ error, tau ]
1190         #   total_iters = iter
1191         #   break
1192          
1193    
1194        if ( rho == 0.0 ):
1195          raise 'TFQMR breakdown, rho=0'
1196        
1197    
1198        rhon = bilinearform(r,w)
1199        beta = rhon / rho;
1200        rho = rhon;
1201        y1 = w + beta * y2;
1202        u1 = Msolve(Aprod(y1))
1203        v = u1 + beta * ( u2 + beta * v )
1204        error = [ error, tau ]
1205        total_iters = iter
1206        
1207        iter = iter + 1
1208    
1209      return x
1210    
1211    
1212    #############################################
1213    
1214    class ArithmeticTuple(object):
1215       """
1216       tuple supporting inplace update x+=y and scaling x=a*y where x,y is an ArithmeticTuple and a is a float.
1217    
1218       example of usage:
1219    
1220       from esys.escript import Data
1221       from numarray import array
1222       a=Data(...)
1223       b=array([1.,4.])
1224       x=ArithmeticTuple(a,b)
1225       y=5.*x
1226    
1227       """
1228       def __init__(self,*args):
1229           """
1230           initialize object with elements args.
1231    
1232           @param args: tuple of object that support implace add (x+=y) and scaling (x=a*y)
1233           """
1234           self.__items=list(args)
1235    
1236       def __len__(self):
1237           """
1238           number of items
1239    
1240           @return: number of items
1241           @rtype: C{int}
1242           """
1243           return len(self.__items)
1244    
1245       def __getitem__(self,index):
1246           """
1247           get an item
1248    
1249           @param index: item to be returned
1250           @type index: C{int}
1251           @return: item with index C{index}
1252           """
1253           return self.__items.__getitem__(index)
1254    
1255       def __mul__(self,other):
1256           """
1257           scaling from the right
1258    
1259           @param other: scaling factor
1260           @type other: C{float}
1261           @return: itemwise self*other
1262           @rtype: L{ArithmeticTuple}
1263           """
1264           out=[]
1265           try:  
1266               l=len(other)
1267               if l!=len(self):
1268                   raise ValueError,"length of of arguments don't match."
1269               for i in range(l): out.append(self[i]*other[i])
1270           except TypeError:
1271           for i in range(len(self)): out.append(self[i]*other)
1272           return ArithmeticTuple(*tuple(out))
1273    
1274       def __rmul__(self,other):
1275           """
1276           scaling from the left
1277    
1278           @param other: scaling factor
1279           @type other: C{float}
1280           @return: itemwise other*self
1281           @rtype: L{ArithmeticTuple}
1282           """
1283           out=[]
1284           try:  
1285               l=len(other)
1286               if l!=len(self):
1287                   raise ValueError,"length of of arguments don't match."
1288               for i in range(l): out.append(other[i]*self[i])
1289           except TypeError:
1290           for i in range(len(self)): out.append(other*self[i])
1291           return ArithmeticTuple(*tuple(out))
1292    
1293       def __div__(self,other):
1294           """
1295           dividing from the right
1296    
1297           @param other: scaling factor
1298           @type other: C{float}
1299           @return: itemwise self/other
1300           @rtype: L{ArithmeticTuple}
1301           """
1302           return self*(1/other)
1303    
1304       def __rdiv__(self,other):
1305           """
1306           dividing from the left
1307    
1308           @param other: scaling factor
1309           @type other: C{float}
1310           @return: itemwise other/self
1311           @rtype: L{ArithmeticTuple}
1312           """
1313           out=[]
1314           try:  
1315               l=len(other)
1316               if l!=len(self):
1317                   raise ValueError,"length of of arguments don't match."
1318               for i in range(l): out.append(other[i]/self[i])
1319           except TypeError:
1320           for i in range(len(self)): out.append(other/self[i])
1321           return ArithmeticTuple(*tuple(out))
1322      
1323       def __iadd__(self,other):
1324           """
1325           in-place add of other to self
1326    
1327           @param other: increment
1328           @type other: C{ArithmeticTuple}
1329           """
1330           if len(self) != len(other):
1331               raise ValueError,"tuple length must match."
1332           for i in range(len(self)):
1333               self.__items[i]+=other[i]
1334           return self
1335    
1336       def __add__(self,other):
1337           """
1338           add other to self
1339    
1340           @param other: increment
1341           @type other: C{ArithmeticTuple}
1342           """
1343           out=[]
1344           try:  
1345               l=len(other)
1346               if l!=len(self):
1347                   raise ValueError,"length of of arguments don't match."
1348               for i in range(l): out.append(self[i]+other[i])
1349           except TypeError:
1350           for i in range(len(self)): out.append(self[i]+other)
1351           return ArithmeticTuple(*tuple(out))
1352    
1353       def __sub__(self,other):
1354           """
1355           subtract other from self
1356    
1357           @param other: increment
1358           @type other: C{ArithmeticTuple}
1359           """
1360           out=[]
1361           try:  
1362               l=len(other)
1363               if l!=len(self):
1364                   raise ValueError,"length of of arguments don't match."
1365               for i in range(l): out.append(self[i]-other[i])
1366           except TypeError:
1367           for i in range(len(self)): out.append(self[i]-other)
1368           return ArithmeticTuple(*tuple(out))
1369      
1370       def __isub__(self,other):
1371           """
1372           subtract other from self
1373    
1374           @param other: increment
1375           @type other: C{ArithmeticTuple}
1376           """
1377           if len(self) != len(other):
1378               raise ValueError,"tuple length must match."
1379           for i in range(len(self)):
1380               self.__items[i]-=other[i]
1381           return self
1382    
1383       def __neg__(self):
1384           """
1385           negate
1386    
1387           """
1388           out=[]
1389           for i in range(len(self)):
1390               out.append(-self[i])
1391           return ArithmeticTuple(*tuple(out))
1392    
1393    
1394    class HomogeneousSaddlePointProblem(object):
1395          """
1396          This provides a framwork for solving linear homogeneous saddle point problem of the form
1397    
1398                 Av+B^*p=f
1399                 Bv    =0
1400    
1401          for the unknowns v and p and given operators A and B and given right hand side f.
1402          B^* is the adjoint operator of B is the given inner product.
1403    
1404          """
1405          def __init__(self,**kwargs):
1406            self.setTolerance()
1407            self.setToleranceReductionFactor()
1408    
1409          def initialize(self):
1410            """
1411            initialize the problem (overwrite)
1412            """
1413            pass
1414          def B(self,v):
1415             """
1416             returns Bv (overwrite)
1417             @rtype: equal to the type of p
1418    
1419             @note: boundary conditions on p should be zero!
1420             """
1421             pass
1422    
1423          def inner(self,p0,p1):
1424             """
1425             returns inner product of two element p0 and p1  (overwrite)
1426            
1427             @type p0: equal to the type of p
1428             @type p1: equal to the type of p
1429             @rtype: equal to the type of p
1430             """
1431             pass
1432    
1433          def solve_A(self,u,p):
1434             """
1435             solves Av=f-Au-B^*p with accuracy self.getReducedTolerance() (overwrite)
1436    
1437             @rtype: equal to the type of v
1438             @note: boundary conditions on v should be zero!
1439             """
1440             pass
1441    
1442          def solve_prec(self,p):
1443             """
1444             provides a preconditioner for BA^{-1}B^* with accuracy self.getReducedTolerance() (overwrite)
1445    
1446             @rtype: equal to the type of p
1447             """
1448             pass
1449    
1450          def stoppingcriterium(self,Bv,v,p):
1451             """
1452             returns a True if iteration is terminated. (overwrite)
1453    
1454             @rtype: C{bool}
1455             """
1456             pass
1457                
1458          def __inner(self,p,r):
1459             return self.inner(p,r[1])
1460    
1461          def __inner_p(self,p1,p2):
1462             return self.inner(p1,p2)
1463          
1464          def __inner_a(self,a1,a2):
1465             return self.inner_a(a1,a2)
1466    
1467          def __inner_a1(self,a1,a2):
1468             return self.inner(a1[1],a2[1])
1469    
1470          def __stoppingcriterium(self,norm_r,r,p):
1471              return self.stoppingcriterium(r[1],r[0],p)
1472    
1473          def __stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None):
1474              return self.stoppingcriterium2(norm_r,norm_b,solver,TOL)
1475    
1476          def setTolerance(self,tolerance=1.e-8):
1477                  self.__tol=tolerance
1478          def getTolerance(self):
1479                  return self.__tol
1480          def setToleranceReductionFactor(self,reduction=0.01):
1481                  self.__reduction=reduction
1482          def getSubProblemTolerance(self):
1483                  return self.__reduction*self.getTolerance()
1484    
1485          def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='PCG',iter_restart=20):
1486                  """
1487                  solves the saddle point problem using initial guesses v and p.
1488    
1489                  @param max_iter: maximum number of iteration steps.
1490                  """
1491                  self.verbose=verbose
1492                  self.show_details=show_details and self.verbose
1493    
1494                  # assume p is known: then v=A^-1(f-B^*p)
1495                  # which leads to BA^-1B^*p = BA^-1f  
1496    
1497              # Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask)      
1498              self.__z=v+self.solve_A(v,p*0)
1499                  Bz=self.B(self.__z)
1500                  #
1501              #   solve BA^-1B^*p = Bz
1502                  #
1503                  #
1504                  #
1505                  self.iter=0
1506              if solver=='GMRES':      
1507                    if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter
1508                    p=GMRES(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.,iter_restart=iter_restart)
1509                    # solve Au=f-B^*p
1510                    #       A(u-v)=f-B^*p-Av
1511                    #       u=v+(u-v)
1512            u=v+self.solve_A(v,p)
1513    
1514              if solver=='TFQMR':      
1515                    if self.verbose: print "enter TFQMR method (iter_max=%s)"%max_iter
1516                    p=TFQMR(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)
1517                    # solve Au=f-B^*p
1518                    #       A(u-v)=f-B^*p-Av
1519                    #       u=v+(u-v)
1520            u=v+self.solve_A(v,p)
1521    
1522              if solver=='MINRES':      
1523                    if self.verbose: print "enter MINRES method (iter_max=%s)"%max_iter
1524                    p=MINRES(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)
1525                    # solve Au=f-B^*p
1526                    #       A(u-v)=f-B^*p-Av
1527                    #       u=v+(u-v)
1528            u=v+self.solve_A(v,p)
1529                  
1530              if solver=='GMRESC':      
1531                    if self.verbose: print "enter GMRES coupled method (iter_max=%s)"%max_iter
1532                    p0=self.solve_prec1(Bz)
1533                #alfa=(1/self.vol)*util.integrate(util.interpolate(p,escript.Function(self.domain)))
1534                    #p-=alfa
1535                    x=GMRES(ArithmeticTuple(self.__z*1.,p0*1),self.__Anext,self.__Mempty,self.__inner_a,self.__stoppingcriterium2,iter_max=max_iter, x=ArithmeticTuple(v*1,p*1),iter_restart=20)
1536                    #x=NewtonGMRES(ArithmeticTuple(self.__z*1.,p0*1),self.__Aprod_Newton2,self.__Mempty,self.__inner_a,self.__stoppingcriterium2,iter_max=max_iter, x=ArithmeticTuple(v*1,p*1),atol=0,rtol=self.getTolerance())
1537    
1538                    # solve Au=f-B^*p
1539                    #       A(u-v)=f-B^*p-Av
1540                    #       u=v+(u-v)
1541                p=x[1]
1542            u=v+self.solve_A(v,p)      
1543            #p=x[1]
1544            #u=x[0]
1545    
1546                  if solver=='PCG':
1547                    #   note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv
1548                    #
1549                    #   with                    Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask)
1550                    #                           A(v-z)= f -Az - B^*p (v-z=0 on fixed_u_mask)
1551                    if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter
1552                    p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p)
1553                u=r[0]  
1554                    # print "DIFF=",util.integrate(p)
1555    
1556                  # print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)
1557    
1558              return u,p
1559    
1560          def __Msolve(self,r):
1561              return self.solve_prec(r[1])
1562    
1563          def __Msolve2(self,r):
1564              return self.solve_prec(r*1)
1565    
1566          def __Mempty(self,r):
1567              return r
1568    
1569    
1570          def __Aprod(self,p):
1571              # return BA^-1B*p
1572              #solve Av =B^*p as Av =f-Az-B^*(-p)
1573              v=self.solve_A(self.__z,-p)
1574              return ArithmeticTuple(v, self.B(v))
1575    
1576          def __Anext(self,x):
1577              # return next v,p
1578              #solve Av =-B^*p as Av =f-Az-B^*p
1579    
1580          pc=x[1]
1581              v=self.solve_A(self.__z,-pc)
1582          p=self.solve_prec1(self.B(v))
1583    
1584              return ArithmeticTuple(v,p)
1585    
1586    
1587          def __Aprod2(self,p):
1588              # return BA^-1B*p
1589              #solve Av =B^*p as Av =f-Az-B^*(-p)
1590          v=self.solve_A(self.__z,-p)
1591              return self.B(v)
1592    
1593          def __Aprod_Newton(self,p):
1594              # return BA^-1B*p - Bz
1595              #solve Av =-B^*p as Av =f-Az-B^*p
1596          v=self.solve_A(self.__z,-p)
1597              return self.B(v-self.__z)
1598    
1599          def __Aprod_Newton2(self,x):
1600              # return BA^-1B*p - Bz
1601              #solve Av =-B^*p as Av =f-Az-B^*p
1602              pc=x[1]
1603          v=self.solve_A(self.__z,-pc)
1604              p=self.solve_prec1(self.B(v-self.__z))
1605              return ArithmeticTuple(v,p)
1606    
1607    
1608    def MaskFromBoundaryTag(domain,*tags):
1609       """
1610       creates a mask on the Solution(domain) function space which one for samples
1611       that touch the boundary tagged by tags.
1612    
1613       usage: m=MaskFromBoundaryTag(domain,"left", "right")
1614    
1615       @param domain: a given domain
1616       @type domain: L{escript.Domain}
1617       @param tags: boundray tags
1618       @type tags: C{str}
1619       @return: a mask which marks samples that are touching the boundary tagged by any of the given tags.
1620       @rtype: L{escript.Data} of rank 0
1621       """
1622       pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1623       d=escript.Scalar(0.,escript.FunctionOnBoundary(domain))
1624       for t in tags: d.setTaggedValue(t,1.)
1625       pde.setValue(y=d)
1626       return util.whereNonZero(pde.getRightHandSide())
1627    #============================================================================================================================================
1628    class SaddlePointProblem(object):
1629       """
1630       This implements a solver for a saddlepoint problem
1631    
1632       M{f(u,p)=0}
1633       M{g(u)=0}
1634    
1635       for u and p. The problem is solved with an inexact Uszawa scheme for p:
1636    
1637       M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})}
1638       M{Q_g (p^{k+1}-p^{k}) =   g(u^{k+1})}
1639    
1640       where Q_f is an approximation of the Jacobiean A_f of f with respect to u  and Q_f is an approximation of
1641       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'
1642       Q_g can be difficult, non-linear conjugate gradient method is applied to solve for p, so Q_g plays
1643       in fact the role of a preconditioner.
1644       """
1645       def __init__(self,verbose=False,*args):
1646           """
1647           initializes the problem
1648    
1649           @param verbose: switches on the printing out some information
1650           @type verbose: C{bool}
1651           @note: this method may be overwritten by a particular saddle point problem
1652           """
1653           print "SaddlePointProblem should not be used anymore!"
1654           if not isinstance(verbose,bool):
1655                raise TypeError("verbose needs to be of type bool.")
1656           self.__verbose=verbose
1657           self.relaxation=1.
1658           DeprecationWarning("SaddlePointProblem should not be used anymore.")
1659    
1660       def trace(self,text):
1661           """
1662           prints text if verbose has been set
1663    
1664           @param text: a text message
1665           @type text: C{str}
1666           """
1667           if self.__verbose: print "%s: %s"%(str(self),text)
1668    
1669       def solve_f(self,u,p,tol=1.e-8):
1670           """
1671           solves
1672    
1673           A_f du = f(u,p)
1674    
1675           with tolerance C{tol} and return du. A_f is Jacobiean of f with respect to u.
1676    
1677           @param u: current approximation of u
1678           @type u: L{escript.Data}
1679           @param p: current approximation of p
1680           @type p: L{escript.Data}
1681           @param tol: tolerance expected for du
1682           @type tol: C{float}
1683           @return: increment du
1684           @rtype: L{escript.Data}
1685           @note: this method has to be overwritten by a particular saddle point problem
1686           """
1687           pass
1688    
1689       def solve_g(self,u,tol=1.e-8):
1690           """
1691           solves
1692    
1693           Q_g dp = g(u)
1694    
1695           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.
1696    
1697           @param u: current approximation of u
1698           @type u: L{escript.Data}
1699           @param tol: tolerance expected for dp
1700           @type tol: C{float}
1701           @return: increment dp
1702           @rtype: L{escript.Data}
1703           @note: this method has to be overwritten by a particular saddle point problem
1704           """
1705           pass
1706    
1707       def inner(self,p0,p1):
1708           """
1709           inner product of p0 and p1 approximating p. Typically this returns integrate(p0*p1)
1710           @return: inner product of p0 and p1
1711           @rtype: C{float}
1712           """
1713           pass
1714    
1715       subiter_max=3
1716       def solve(self,u0,p0,tolerance=1.e-6,tolerance_u=None,iter_max=100,accepted_reduction=0.995,relaxation=None):
1717            """
1718            runs the solver
1719    
1720            @param u0: initial guess for C{u}
1721            @type u0: L{esys.escript.Data}
1722            @param p0: initial guess for C{p}
1723            @type p0: L{esys.escript.Data}
1724            @param tolerance: tolerance for relative error in C{u} and C{p}
1725            @type tolerance: positive C{float}
1726            @param tolerance_u: tolerance for relative error in C{u} if different from C{tolerance}
1727            @type tolerance_u: positive C{float}
1728            @param iter_max: maximum number of iteration steps.
1729            @type iter_max: C{int}
1730            @param accepted_reduction: if the norm  g cannot be reduced by C{accepted_reduction} backtracking to adapt the
1731                                       relaxation factor. If C{accepted_reduction=None} no backtracking is used.
1732            @type accepted_reduction: positive C{float} or C{None}
1733            @param relaxation: initial relaxation factor. If C{relaxation==None}, the last relaxation factor is used.
1734            @type relaxation: C{float} or C{None}
1735            """
1736            tol=1.e-2
1737            if tolerance_u==None: tolerance_u=tolerance
1738            if not relaxation==None: self.relaxation=relaxation
1739            if accepted_reduction ==None:
1740                  angle_limit=0.
1741            elif accepted_reduction>=1.:
1742                  angle_limit=0.
1743            else:
1744                  angle_limit=util.sqrt(1-accepted_reduction**2)
1745            self.iter=0
1746            u=u0
1747            p=p0
1748            #
1749            #   initialize things:
1750            #
1751            converged=False
1752            #
1753            #  start loop:
1754            #
1755            #  initial search direction is g
1756            #
1757            while not converged :
1758                if self.iter>iter_max:
1759                    raise ArithmeticError("no convergence after %s steps."%self.iter)
1760                f_new=self.solve_f(u,p,tol)
1761                norm_f_new = util.Lsup(f_new)
1762                u_new=u-f_new
1763                g_new=self.solve_g(u_new,tol)
1764                self.iter+=1
1765                norm_g_new = util.sqrt(self.inner(g_new,g_new))
1766                if norm_f_new==0. and norm_g_new==0.: return u, p
1767                if self.iter>1 and not accepted_reduction==None:
1768                   #
1769                   #   did we manage to reduce the norm of G? I
1770                   #   if not we start a backtracking procedure
1771                   #
1772                   # print "new/old norm = ",norm_g_new, norm_g, norm_g_new/norm_g
1773                   if norm_g_new > accepted_reduction * norm_g:
1774                      sub_iter=0
1775                      s=self.relaxation
1776                      d=g
1777                      g_last=g
1778                      self.trace("    start substepping: f = %s, g = %s, relaxation = %s."%(norm_f_new, norm_g_new, s))
1779                      while sub_iter < self.subiter_max and  norm_g_new > accepted_reduction * norm_g:
1780                         dg= g_new-g_last
1781                         norm_dg=abs(util.sqrt(self.inner(dg,dg))/self.relaxation)
1782                         rad=self.inner(g_new,dg)/self.relaxation
1783                         # print "   ",sub_iter,": rad, norm_dg:",abs(rad), norm_dg*norm_g_new * angle_limit
1784                         # print "   ",sub_iter,": rad, norm_dg:",rad, norm_dg, norm_g_new, norm_g
1785                         if abs(rad) < norm_dg*norm_g_new * angle_limit:
1786                             if sub_iter>0: self.trace("    no further improvements expected from backtracking.")
1787                             break
1788                         r=self.relaxation
1789                         self.relaxation= - rad/norm_dg**2
1790                         s+=self.relaxation
1791                         #####
1792                         # a=g_new+self.relaxation*dg/r
1793                         # print "predicted new norm = ",util.sqrt(self.inner(a,a)),util.sqrt(self.inner(g_new,g_new)), self.relaxation
1794                         #####
1795                         g_last=g_new
1796                         p+=self.relaxation*d
1797                         f_new=self.solve_f(u,p,tol)
1798                         u_new=u-f_new
1799                         g_new=self.solve_g(u_new,tol)
1800                         self.iter+=1
1801                         norm_f_new = util.Lsup(f_new)
1802                         norm_g_new = util.sqrt(self.inner(g_new,g_new))
1803                         # print "   ",sub_iter," new g norm",norm_g_new
1804                         self.trace("    %s th sub-step: f = %s, g = %s, relaxation = %s."%(sub_iter, norm_f_new, norm_g_new, s))
1805                         #
1806                         #   can we expect reduction of g?
1807                         #
1808                         # u_last=u_new
1809                         sub_iter+=1
1810                      self.relaxation=s
1811                #
1812                #  check for convergence:
1813                #
1814                norm_u_new = util.Lsup(u_new)
1815                p_new=p+self.relaxation*g_new
1816                norm_p_new = util.sqrt(self.inner(p_new,p_new))
1817                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))
1818    
1819                if self.iter>1:
1820                   dg2=g_new-g
1821                   df2=f_new-f
1822                   norm_dg2=util.sqrt(self.inner(dg2,dg2))
1823                   norm_df2=util.Lsup(df2)
1824                   # print norm_g_new, norm_g, norm_dg, norm_p, tolerance
1825                   tol_eq_g=tolerance*norm_dg2/(norm_g*abs(self.relaxation))*norm_p_new
1826                   tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new
1827                   if norm_g_new <= tol_eq_g and norm_f_new <= tol_eq_f:
1828                       converged=True
1829                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
1830            self.trace("convergence after %s steps."%self.iter)
1831            return u,p

Legend:
Removed from v.351  
changed lines
  Added in v.2061

  ViewVC Help
Powered by ViewVC 1.1.26