/[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 155 by jgs, Wed Nov 9 02:02:19 2005 UTC revision 1956 by gross, Mon Nov 3 05:08:42 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.
24    
25  Currently includes:  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
28        - 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:
54      """
55      a simple mechanism to manage time dependend values.
56    
57      typical usage is::
58    
59         dt=0.1 # time increment
60         tm=TimeIntegrationManager(inital_value,p=1)
61         while t<1.
62             v_guess=tm.extrapolate(dt) # extrapolate to t+dt
63             v=...
64             tm.checkin(dt,v)
65             t+=dt
66    
67      @note: currently only p=1 is supported.
68      """
69      def __init__(self,*inital_values,**kwargs):
70         """
71         sets up the value manager where inital_value is the initial value and p is order used for extrapolation
72         """
73         if kwargs.has_key("p"):
74                self.__p=kwargs["p"]
75         else:
76                self.__p=1
77         if kwargs.has_key("time"):
78                self.__t=kwargs["time"]
79         else:
80                self.__t=0.
81         self.__v_mem=[inital_values]
82         self.__order=0
83         self.__dt_mem=[]
84         self.__num_val=len(inital_values)
85    
86      def getTime(self):
87          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):
96          """
97          adds new values to the manager. the p+1 last value get lost
98          """
99          o=min(self.__order+1,self.__p)
100          self.__order=min(self.__order+1,self.__p)
101          v_mem_new=[values]
102          dt_mem_new=[dt]
103          for i in range(o-1):
104             v_mem_new.append(self.__v_mem[i])
105             dt_mem_new.append(self.__dt_mem[i])
106          v_mem_new.append(self.__v_mem[o-1])
107          self.__order=o
108          self.__v_mem=v_mem_new
109          self.__dt_mem=dt_mem_new
110          self.__t+=dt
111    
112      def extrapolate(self,dt):
113          """
114          extrapolates to dt forward in time.
115          """
116          if self.__order==0:
117             out=self.__v_mem[0]
118          else:
119            out=[]
120            for i in range(self.__num_val):
121               out.append((1.+dt/self.__dt_mem[0])*self.__v_mem[0][i]-dt/self.__dt_mem[0]*self.__v_mem[1][i])
122    
123          if len(out)==0:
124             return None
125          elif len(out)==1:
126             return out[0]
127          else:
128             return out
129    
130    
131  class Projector:  class Projector:
132    """    """
# Line 33  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 70  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 86  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 105  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 131  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}. The returned object needs to be of the same type like argument C{b}.
517       @param Msolve: solves Mx=r
518       @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
519    type like argument C{x}.
520       @param bilinearform: inner product C{<x,r>}
521       @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}.
522       @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.
523       @type stoppingcriterium: function that returns C{True} or C{False}
524       @param x: an initial guess for the solution. If no C{x} is given 0*b is used.
525       @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
526       @param iter_max: maximum number of iteration steps.
527       @type iter_max: C{int}
528       @return: the solution approximation and the corresponding residual
529       @rtype: C{tuple}
530       @warning: C{b} and C{x} are altered.
531       """
532       iter=0
533       if x==None:
534          x=0*b
535       else:
536          b += (-1)*Aprod(x)
537       r=b
538       rhat=Msolve(r)
539       d = rhat
540       rhat_dot_r = bilinearform(rhat, r)
541       if rhat_dot_r<0: raise NegativeNorm,"negative norm."
542    
543       while not stoppingcriterium(math.sqrt(rhat_dot_r),r,x):
544           iter+=1
545           if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
546    
547           q=Aprod(d)
548           alpha = rhat_dot_r / bilinearform(d, q)
549           x += alpha * d
550           r += (-alpha) * q
551    
552           rhat=Msolve(r)
553           rhat_dot_r_new = bilinearform(rhat, r)
554           beta = rhat_dot_r_new / rhat_dot_r
555           rhat+=beta * d
556           d=rhat
557    
558           rhat_dot_r = rhat_dot_r_new
559           if rhat_dot_r<0: raise NegativeNorm,"negative norm."
560    
561       return x,r
562    
563    class Defect(object):
564        """
565        defines a non-linear defect F(x) of a variable x
566        """
567        def __init__(self):
568            """
569            initialize defect
570            """
571            self.setDerivativeIncrementLength()
572    
573        def bilinearform(self, x0, x1):
574            """
575            returns the inner product of x0 and x1
576            @param x0: a value for x
577            @param x1: a value for x
578            @return: the inner product of x0 and x1
579            @rtype: C{float}
580            """
581            return 0
582          
583        def norm(self,x):
584            """
585            the norm of argument C{x}
586    
587            @param x: a value for x
588            @return: norm of argument x
589            @rtype: C{float}
590            @note: by default C{sqrt(self.bilinearform(x,x)} is retrurned.
591            """
592            s=self.bilinearform(x,x)
593            if s<0: raise NegativeNorm,"negative norm."
594            return math.sqrt(s)
595    
596    
597        def eval(self,x):
598            """
599            returns the value F of a given x
600    
601            @param x: value for which the defect C{F} is evalulated.
602            @return: value of the defect at C{x}
603            """
604            return 0
605    
606        def __call__(self,x):
607            return self.eval(x)
608    
609        def setDerivativeIncrementLength(self,inc=math.sqrt(util.EPSILON)):
610            """
611            sets the relative length of the increment used to approximate the derivative of the defect
612            the increment is inc*norm(x)/norm(v)*v in the direction of v with x as a staring point.
613    
614            @param inc: relative increment length
615            @type inc: positive C{float}
616            """
617            if inc<=0: raise ValueError,"positive increment required."
618            self.__inc=inc
619    
620        def getDerivativeIncrementLength(self):
621            """
622            returns the relative increment length used to approximate the derivative of the defect
623            @return: value of the defect at C{x}
624            @rtype: positive C{float}
625            """
626            return self.__inc
627    
628        def derivative(self, F0, x0, v, v_is_normalised=True):
629            """
630            returns the directional derivative at x0 in the direction of v
631    
632            @param F0: value of this defect at x0
633            @param x0: value at which derivative is calculated.
634            @param v: direction
635            @param v_is_normalised: is true to indicate that C{v} is nomalized (self.norm(v)=0)
636            @return: derivative of this defect at x0 in the direction of C{v}
637            @note: by default numerical evaluation (self.eval(x0+eps*v)-F0)/eps is used but this method
638            maybe oepsnew verwritten to use exact evalution.
639            """
640            normx=self.norm(x0)
641            if normx>0:
642                 epsnew = self.getDerivativeIncrementLength() * normx
643            else:
644                 epsnew = self.getDerivativeIncrementLength()
645            if not v_is_normalised:
646                normv=self.norm(v)
647                if normv<=0:
648                   return F0*0
649                else:
650                   epsnew /= normv
651            F1=self.eval(x0 + epsnew * v)
652            return (F1-F0)/epsnew
653    
654    ######################################    
655    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):
656       """
657       solves a non-linear problem M{F(x)=0} for unknown M{x} using the stopping criterion:
658    
659       M{norm(F(x) <= atol + rtol * norm(F(x0)}
660      
661       where M{x0} is the initial guess.
662    
663       @param defect: object defining the the function M{F}, C{defect.norm} defines the M{norm} used in the stopping criterion.
664       @type defect: L{Defect}
665       @param x: initial guess for the solution, C{x} is altered.
666       @type x: any object type allowing basic operations such as  L{numarray.NumArray}, L{Data}
667       @param iter_max: maximum number of iteration steps
668       @type iter_max: positive C{int}
669       @param sub_iter_max:
670       @type sub_iter_max:
671       @param atol: absolute tolerance for the solution
672       @type atol: positive C{float}
673       @param rtol: relative tolerance for the solution
674       @type rtol: positive C{float}
675       @param gamma: tolerance safety factor for inner iteration
676       @type gamma: positive C{float}, less than 1
677       @param sub_tol_max: upper bound for inner tolerance.
678       @type sub_tol_max: positive C{float}, less than 1
679       @return: an approximation of the solution with the desired accuracy
680       @rtype: same type as the initial guess.
681       """
682       lmaxit=iter_max
683       if atol<0: raise ValueError,"atol needs to be non-negative."
684       if rtol<0: raise ValueError,"rtol needs to be non-negative."
685       if rtol+atol<=0: raise ValueError,"rtol or atol needs to be non-negative."
686       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
687       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
688    
689       F=defect(x)
690       fnrm=defect.norm(F)
691       stop_tol=atol + rtol*fnrm
692       sub_tol=sub_tol_max
693       if verbose: print "NewtonGMRES: initial residual = %e."%fnrm
694       if verbose: print "             tolerance = %e."%sub_tol
695       iter=1
696       #
697       # main iteration loop
698       #
699       while not fnrm<=stop_tol:
700                if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
701                #
702            #   adjust sub_tol_
703            #
704                if iter > 1:
705               rat=fnrm/fnrmo
706                   sub_tol_old=sub_tol
707               sub_tol=gamma*rat**2
708               if gamma*sub_tol_old**2 > .1: sub_tol=max(sub_tol,gamma*sub_tol_old**2)
709               sub_tol=max(min(sub_tol,sub_tol_max), .5*stop_tol/fnrm)
710            #
711            # calculate newton increment xc
712                #     if iter_max in __FDGMRES is reached MaxIterReached is thrown
713                #     if iter_restart -1 is returned as sub_iter
714                #     if  atol is reached sub_iter returns the numer of steps performed to get there
715                #
716                #  
717                if verbose: print "             subiteration (GMRES) is called with relative tolerance %e."%sub_tol
718                try:
719                   xc, sub_iter=__FDGMRES(F, defect, x, sub_tol*fnrm, iter_max=iter_max-iter, iter_restart=sub_iter_max)
720                except MaxIterReached:
721                   raise MaxIterReached,"maximum number of %s steps reached."%iter_max
722                if sub_iter<0:
723                   iter+=sub_iter_max
724                else:
725                   iter+=sub_iter
726                # ====
727            x+=xc
728                F=defect(x)
729            iter+=1
730                fnrmo, fnrm=fnrm, defect.norm(F)
731                if verbose: print "             step %s: residual %e."%(iter,fnrm)
732       if verbose: print "NewtonGMRES: completed after %s steps."%iter
733       return x
734    
735    def __givapp(c,s,vin):
736        """
737        apply a sequence of Givens rotations (c,s) to the recuirsively to the vector vin
738        @warning: C{vin} is altered.
739        """
740        vrot=vin
741        if isinstance(c,float):
742            vrot=[c*vrot[0]-s*vrot[1],s*vrot[0]+c*vrot[1]]
743        else:
744            for i in range(len(c)):
745                w1=c[i]*vrot[i]-s[i]*vrot[i+1]
746            w2=s[i]*vrot[i]+c[i]*vrot[i+1]
747                vrot[i:i+2]=w1,w2
748        return vrot
749    
750    def __FDGMRES(F0, defect, x0, atol, iter_max=100, iter_restart=20):
751       h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)
752       c=numarray.zeros(iter_restart,numarray.Float64)
753       s=numarray.zeros(iter_restart,numarray.Float64)
754       g=numarray.zeros(iter_restart,numarray.Float64)
755       v=[]
756    
757       rho=defect.norm(F0)
758       if rho<=0.: return x0*0
759      
760       v.append(-F0/rho)
761       g[0]=rho
762       iter=0
763       while rho > atol and iter<iter_restart-1:
764    
765        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
766    
767            p=defect.derivative(F0,x0,v[iter], v_is_normalised=True)
768        v.append(p)
769    
770        v_norm1=defect.norm(v[iter+1])
771    
772            # Modified Gram-Schmidt
773        for j in range(iter+1):
774             h[j][iter]=defect.bilinearform(v[j],v[iter+1])  
775             v[iter+1]-=h[j][iter]*v[j]
776          
777        h[iter+1][iter]=defect.norm(v[iter+1])
778        v_norm2=h[iter+1][iter]
779    
780            # Reorthogonalize if needed
781        if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
782            for j in range(iter+1):  
783               hr=defect.bilinearform(v[j],v[iter+1])
784                   h[j][iter]=h[j][iter]+hr
785                   v[iter+1] -= hr*v[j]
786    
787            v_norm2=defect.norm(v[iter+1])
788            h[iter+1][iter]=v_norm2
789            #   watch out for happy breakdown
790            if not v_norm2 == 0:
791                    v[iter+1]=v[iter+1]/h[iter+1][iter]
792    
793            #   Form and store the information for the new Givens rotation
794        if iter > 0 :
795            hhat=numarray.zeros(iter+1,numarray.Float64)
796            for i in range(iter+1) : hhat[i]=h[i][iter]
797            hhat=__givapp(c[0:iter],s[0:iter],hhat);
798                for i in range(iter+1) : h[i][iter]=hhat[i]
799    
800        mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])
801    
802        if mu!=0 :
803            c[iter]=h[iter][iter]/mu
804            s[iter]=-h[iter+1][iter]/mu
805            h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter]
806            h[iter+1][iter]=0.0
807            g[iter:iter+2]=__givapp(c[iter],s[iter],g[iter:iter+2])
808    
809            # Update the residual norm
810            rho=abs(g[iter+1])
811        iter+=1
812    
813       # At this point either iter > iter_max or rho < tol.
814       # It's time to compute x and leave.        
815       if iter > 0 :
816         y=numarray.zeros(iter,numarray.Float64)    
817         y[iter-1] = g[iter-1] / h[iter-1][iter-1]
818         if iter > 1 :  
819            i=iter-2  
820            while i>=0 :
821              y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]
822              i=i-1
823         xhat=v[iter-1]*y[iter-1]
824         for i in range(iter-1):
825        xhat += v[i]*y[i]
826       else :
827          xhat=v[0] * 0
828    
829       if iter<iter_restart-1:
830          stopped=iter
831       else:
832          stopped=-1
833    
834       return xhat,stopped
835    
836    ##############################################
837    def GMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):
838    ################################################
839       m=iter_restart
840       iter=0
841       xc=x
842       while True:
843          if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached"%iter_max
844          xc,stopped=__GMRESm(b*1, Aprod, Msolve, bilinearform, stoppingcriterium, x=xc*1, iter_max=iter_max-iter, iter_restart=m)
845          if stopped: break
846          iter+=iter_restart    
847       return xc
848    
849    def __GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):
850       iter=0
851       r=Msolve(b)
852       r_dot_r = bilinearform(r, r)
853       if r_dot_r<0: raise NegativeNorm,"negative norm."
854       norm_b=math.sqrt(r_dot_r)
855    
856       if x==None:
857          x=0*b
858       else:
859          r=Msolve(b-Aprod(x))
860          r_dot_r = bilinearform(r, r)
861          if r_dot_r<0: raise NegativeNorm,"negative norm."
862      
863       h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)
864       c=numarray.zeros(iter_restart,numarray.Float64)
865       s=numarray.zeros(iter_restart,numarray.Float64)
866       g=numarray.zeros(iter_restart,numarray.Float64)
867       v=[]
868    
869       rho=math.sqrt(r_dot_r)
870      
871       v.append(r/rho)
872       g[0]=rho
873    
874       while not (stoppingcriterium(rho,norm_b) or iter==iter_restart-1):
875    
876        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
877    
878        p=Msolve(Aprod(v[iter]))
879    
880        v.append(p)
881    
882        v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
883    
884    # Modified Gram-Schmidt
885        for j in range(iter+1):
886          h[j][iter]=bilinearform(v[j],v[iter+1])  
887          v[iter+1]-=h[j][iter]*v[j]
888          
889        h[iter+1][iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
890        v_norm2=h[iter+1][iter]
891    
892    # Reorthogonalize if needed
893        if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
894         for j in range(iter+1):  
895            hr=bilinearform(v[j],v[iter+1])
896                h[j][iter]=h[j][iter]+hr
897                v[iter+1] -= hr*v[j]
898    
899         v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
900         h[iter+1][iter]=v_norm2
901    
902    #   watch out for happy breakdown
903            if not v_norm2 == 0:
904             v[iter+1]=v[iter+1]/h[iter+1][iter]
905    
906    #   Form and store the information for the new Givens rotation
907        if iter > 0 :
908            hhat=numarray.zeros(iter+1,numarray.Float64)
909            for i in range(iter+1) : hhat[i]=h[i][iter]
910            hhat=__givapp(c[0:iter],s[0:iter],hhat);
911                for i in range(iter+1) : h[i][iter]=hhat[i]
912    
913        mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])
914    
915        if mu!=0 :
916            c[iter]=h[iter][iter]/mu
917            s[iter]=-h[iter+1][iter]/mu
918            h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter]
919            h[iter+1][iter]=0.0
920            g[iter:iter+2]=__givapp(c[iter],s[iter],g[iter:iter+2])
921    
922    # Update the residual norm
923                  
924            rho=abs(g[iter+1])
925        iter+=1
926    
927    # At this point either iter > iter_max or rho < tol.
928    # It's time to compute x and leave.        
929    
930       if iter > 0 :
931         y=numarray.zeros(iter,numarray.Float64)    
932         y[iter-1] = g[iter-1] / h[iter-1][iter-1]
933         if iter > 1 :  
934            i=iter-2  
935            while i>=0 :
936              y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]
937              i=i-1
938         xhat=v[iter-1]*y[iter-1]
939         for i in range(iter-1):
940        xhat += v[i]*y[i]
941       else : xhat=v[0]
942    
943       x += xhat
944       if iter<iter_restart-1:
945          stopped=True
946       else:
947          stopped=False
948    
949       return x,stopped
950    
951    #################################################
952    def MINRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
953    #################################################
954        #
955        #  minres solves the system of linear equations Ax = b
956        #  where A is a symmetric matrix (possibly indefinite or singular)
957        #  and b is a given vector.
958        #  
959        #  "A" may be a dense or sparse matrix (preferably sparse!)
960        #  or the name of a function such that
961        #               y = A(x)
962        #  returns the product y = Ax for any given vector x.
963        #
964        #  "M" defines a positive-definite preconditioner M = C C'.
965        #  "M" may be a dense or sparse matrix (preferably sparse!)
966        #  or the name of a function such that
967        #  solves the system My = x for any given vector x.
968        #
969        #
970        
971        #------------------------------------------------------------------
972        # Set up y and v for the first Lanczos vector v1.
973        # y  =  beta1 P' v1,  where  P = C**(-1).
974        # v is really P' v1.
975        #------------------------------------------------------------------
976        if x==None:
977          x=0*b
978        else:
979          b += (-1)*Aprod(x)
980    
981        r1    = b
982        y = Msolve(b)
983        beta1 = bilinearform(y,b)
984    
985        if beta1< 0: raise NegativeNorm,"negative norm."
986    
987        #  If b = 0 exactly, stop with x = 0.
988        if beta1==0: return x*0.
989    
990        if beta1> 0:
991          beta1  = math.sqrt(beta1)      
992    
993        #------------------------------------------------------------------
994        # Initialize quantities.
995        # ------------------------------------------------------------------
996        iter   = 0
997        Anorm = 0
998        ynorm = 0
999        oldb   = 0
1000        beta   = beta1
1001        dbar   = 0
1002        epsln  = 0
1003        phibar = beta1
1004        rhs1   = beta1
1005        rhs2   = 0
1006        rnorm  = phibar
1007        tnorm2 = 0
1008        ynorm2 = 0
1009        cs     = -1
1010        sn     = 0
1011        w      = b*0.
1012        w2     = b*0.
1013        r2     = r1
1014        eps    = 0.0001
1015    
1016        #---------------------------------------------------------------------
1017        # Main iteration loop.
1018        # --------------------------------------------------------------------
1019        while not stoppingcriterium(rnorm,Anorm*ynorm,'MINRES'):    #  checks ||r|| < (||A|| ||x||) * TOL
1020    
1021        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
1022            iter    = iter  +  1
1023    
1024            #-----------------------------------------------------------------
1025            # Obtain quantities for the next Lanczos vector vk+1, k = 1, 2,...
1026            # The general iteration is similar to the case k = 1 with v0 = 0:
1027            #
1028            #   p1      = Operator * v1  -  beta1 * v0,
1029            #   alpha1  = v1'p1,
1030            #   q2      = p2  -  alpha1 * v1,
1031            #   beta2^2 = q2'q2,
1032            #   v2      = (1/beta2) q2.
1033            #
1034            # Again, y = betak P vk,  where  P = C**(-1).
1035            #-----------------------------------------------------------------
1036            s = 1/beta                 # Normalize previous vector (in y).
1037            v = s*y                    # v = vk if P = I
1038        
1039            y      = Aprod(v)
1040        
1041            if iter >= 2:
1042              y = y - (beta/oldb)*r1
1043    
1044            alfa   = bilinearform(v,y)              # alphak
1045            y      += (- alfa/beta)*r2
1046            r1     = r2
1047            r2     = y
1048            y = Msolve(r2)
1049            oldb   = beta                           # oldb = betak
1050            beta   = bilinearform(y,r2)             # beta = betak+1^2
1051            if beta < 0: raise NegativeNorm,"negative norm."
1052    
1053            beta   = math.sqrt( beta )
1054            tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta
1055            
1056            if iter==1:                 # Initialize a few things.
1057              gmax   = abs( alfa )      # alpha1
1058              gmin   = gmax             # alpha1
1059    
1060            # Apply previous rotation Qk-1 to get
1061            #   [deltak epslnk+1] = [cs  sn][dbark    0   ]
1062            #   [gbar k dbar k+1]   [sn -cs][alfak betak+1].
1063        
1064            oldeps = epsln
1065            delta  = cs * dbar  +  sn * alfa  # delta1 = 0         deltak
1066            gbar   = sn * dbar  -  cs * alfa  # gbar 1 = alfa1     gbar k
1067            epsln  =               sn * beta  # epsln2 = 0         epslnk+1
1068            dbar   =            -  cs * beta  # dbar 2 = beta2     dbar k+1
1069    
1070            # Compute the next plane rotation Qk
1071    
1072            gamma  = math.sqrt(gbar*gbar+beta*beta)  # gammak
1073            gamma  = max(gamma,eps)
1074            cs     = gbar / gamma             # ck
1075            sn     = beta / gamma             # sk
1076            phi    = cs * phibar              # phik
1077            phibar = sn * phibar              # phibark+1
1078    
1079            # Update  x.
1080    
1081            denom = 1/gamma
1082            w1    = w2
1083            w2    = w
1084            w     = (v - oldeps*w1 - delta*w2) * denom
1085            x     +=  phi*w
1086    
1087            # Go round again.
1088    
1089            gmax   = max(gmax,gamma)
1090            gmin   = min(gmin,gamma)
1091            z      = rhs1 / gamma
1092            ynorm2 = z*z  +  ynorm2
1093            rhs1   = rhs2 -  delta*z
1094            rhs2   =      -  epsln*z
1095    
1096            # Estimate various norms and test for convergence.
1097    
1098            Anorm  = math.sqrt( tnorm2 )
1099            ynorm  = math.sqrt( ynorm2 )
1100    
1101            rnorm  = phibar
1102    
1103        return x
1104    
1105    def TFQMR(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
1106    
1107    # TFQMR solver for linear systems
1108    #
1109    #
1110    # initialization
1111    #
1112      errtol = math.sqrt(bilinearform(b,b))
1113      norm_b=errtol
1114      kmax  = iter_max
1115      error = []
1116    
1117      if math.sqrt(bilinearform(x,x)) != 0.0:
1118        r = b - Aprod(x)
1119      else:
1120        r = b
1121    
1122      r=Msolve(r)
1123    
1124      u1=0
1125      u2=0
1126      y1=0
1127      y2=0
1128    
1129      w = r
1130      y1 = r
1131      iter = 0
1132      d = 0
1133      
1134      v = Msolve(Aprod(y1))
1135      u1 = v
1136      
1137      theta = 0.0;
1138      eta = 0.0;
1139      tau = math.sqrt(bilinearform(r,r))
1140      error = [ error, tau ]
1141      rho = tau * tau
1142      m=1
1143    #
1144    #  TFQMR iteration
1145    #
1146    #  while ( iter < kmax-1 ):
1147      
1148      while not stoppingcriterium(tau*math.sqrt ( m + 1 ),norm_b,'TFQMR'):
1149        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
1150    
1151        sigma = bilinearform(r,v)
1152    
1153        if ( sigma == 0.0 ):
1154          raise 'TFQMR breakdown, sigma=0'
1155        
1156    
1157        alpha = rho / sigma
1158    
1159        for j in range(2):
1160    #
1161    #   Compute y2 and u2 only if you have to
1162    #
1163          if ( j == 1 ):
1164            y2 = y1 - alpha * v
1165            u2 = Msolve(Aprod(y2))
1166    
1167          m = 2 * (iter+1) - 2 + (j+1)
1168          if j==0:
1169             w = w - alpha * u1
1170             d = y1 + ( theta * theta * eta / alpha ) * d
1171          if j==1:
1172             w = w - alpha * u2
1173             d = y2 + ( theta * theta * eta / alpha ) * d
1174    
1175          theta = math.sqrt(bilinearform(w,w))/ tau
1176          c = 1.0 / math.sqrt ( 1.0 + theta * theta )
1177          tau = tau * theta * c
1178          eta = c * c * alpha
1179          x = x + eta * d
1180    #
1181    #  Try to terminate the iteration at each pass through the loop
1182    #
1183         # if ( tau * math.sqrt ( m + 1 ) <= errtol ):
1184         #   error = [ error, tau ]
1185         #   total_iters = iter
1186         #   break
1187          
1188    
1189        if ( rho == 0.0 ):
1190          raise 'TFQMR breakdown, rho=0'
1191        
1192    
1193        rhon = bilinearform(r,w)
1194        beta = rhon / rho;
1195        rho = rhon;
1196        y1 = w + beta * y2;
1197        u1 = Msolve(Aprod(y1))
1198        v = u1 + beta * ( u2 + beta * v )
1199        error = [ error, tau ]
1200        total_iters = iter
1201        
1202        iter = iter + 1
1203    
1204      return x
1205    
1206    
1207    #############################################
1208    
1209    class ArithmeticTuple(object):
1210       """
1211       tuple supporting inplace update x+=y and scaling x=a*y where x,y is an ArithmeticTuple and a is a float.
1212    
1213       example of usage:
1214    
1215       from esys.escript import Data
1216       from numarray import array
1217       a=Data(...)
1218       b=array([1.,4.])
1219       x=ArithmeticTuple(a,b)
1220       y=5.*x
1221    
1222       """
1223       def __init__(self,*args):
1224           """
1225           initialize object with elements args.
1226    
1227           @param args: tuple of object that support implace add (x+=y) and scaling (x=a*y)
1228           """
1229           self.__items=list(args)
1230    
1231       def __len__(self):
1232           """
1233           number of items
1234    
1235           @return: number of items
1236           @rtype: C{int}
1237           """
1238           return len(self.__items)
1239    
1240       def __getitem__(self,index):
1241           """
1242           get an item
1243    
1244           @param index: item to be returned
1245           @type index: C{int}
1246           @return: item with index C{index}
1247           """
1248           return self.__items.__getitem__(index)
1249    
1250       def __mul__(self,other):
1251           """
1252           scaling from the right
1253    
1254           @param other: scaling factor
1255           @type other: C{float}
1256           @return: itemwise self*other
1257           @rtype: L{ArithmeticTuple}
1258           """
1259           out=[]
1260           try:  
1261               l=len(other)
1262               if l!=len(self):
1263                   raise ValueError,"length of of arguments don't match."
1264               for i in range(l): out.append(self[i]*other[i])
1265           except TypeError:
1266           for i in range(len(self)): out.append(self[i]*other)
1267           return ArithmeticTuple(*tuple(out))
1268    
1269       def __rmul__(self,other):
1270           """
1271           scaling from the left
1272    
1273           @param other: scaling factor
1274           @type other: C{float}
1275           @return: itemwise other*self
1276           @rtype: L{ArithmeticTuple}
1277           """
1278           out=[]
1279           try:  
1280               l=len(other)
1281               if l!=len(self):
1282                   raise ValueError,"length of of arguments don't match."
1283               for i in range(l): out.append(other[i]*self[i])
1284           except TypeError:
1285           for i in range(len(self)): out.append(other*self[i])
1286           return ArithmeticTuple(*tuple(out))
1287    
1288       def __div__(self,other):
1289           """
1290           dividing from the right
1291    
1292           @param other: scaling factor
1293           @type other: C{float}
1294           @return: itemwise self/other
1295           @rtype: L{ArithmeticTuple}
1296           """
1297           return self*(1/other)
1298    
1299       def __rdiv__(self,other):
1300           """
1301           dividing from the left
1302    
1303           @param other: scaling factor
1304           @type other: C{float}
1305           @return: itemwise other/self
1306           @rtype: L{ArithmeticTuple}
1307           """
1308           out=[]
1309           try:  
1310               l=len(other)
1311               if l!=len(self):
1312                   raise ValueError,"length of of arguments don't match."
1313               for i in range(l): out.append(other[i]/self[i])
1314           except TypeError:
1315           for i in range(len(self)): out.append(other/self[i])
1316           return ArithmeticTuple(*tuple(out))
1317      
1318       def __iadd__(self,other):
1319           """
1320           in-place add of other to self
1321    
1322           @param other: increment
1323           @type other: C{ArithmeticTuple}
1324           """
1325           if len(self) != len(other):
1326               raise ValueError,"tuple length must match."
1327           for i in range(len(self)):
1328               self.__items[i]+=other[i]
1329           return self
1330    
1331       def __add__(self,other):
1332           """
1333           add other to self
1334    
1335           @param other: increment
1336           @type other: C{ArithmeticTuple}
1337           """
1338           out=[]
1339           try:  
1340               l=len(other)
1341               if l!=len(self):
1342                   raise ValueError,"length of of arguments don't match."
1343               for i in range(l): out.append(self[i]+other[i])
1344           except TypeError:
1345           for i in range(len(self)): out.append(self[i]+other)
1346           return ArithmeticTuple(*tuple(out))
1347    
1348       def __sub__(self,other):
1349           """
1350           subtract other from self
1351    
1352           @param other: increment
1353           @type other: C{ArithmeticTuple}
1354           """
1355           out=[]
1356           try:  
1357               l=len(other)
1358               if l!=len(self):
1359                   raise ValueError,"length of of arguments don't match."
1360               for i in range(l): out.append(self[i]-other[i])
1361           except TypeError:
1362           for i in range(len(self)): out.append(self[i]-other)
1363           return ArithmeticTuple(*tuple(out))
1364      
1365       def __isub__(self,other):
1366           """
1367           subtract other from self
1368    
1369           @param other: increment
1370           @type other: C{ArithmeticTuple}
1371           """
1372           if len(self) != len(other):
1373               raise ValueError,"tuple length must match."
1374           for i in range(len(self)):
1375               self.__items[i]-=other[i]
1376           return self
1377    
1378       def __neg__(self):
1379           """
1380           negate
1381    
1382           """
1383           out=[]
1384           for i in range(len(self)):
1385               out.append(-self[i])
1386           return ArithmeticTuple(*tuple(out))
1387    
1388    
1389    class HomogeneousSaddlePointProblem(object):
1390          """
1391          This provides a framwork for solving linear homogeneous saddle point problem of the form
1392    
1393                 Av+B^*p=f
1394                 Bv    =0
1395    
1396          for the unknowns v and p and given operators A and B and given right hand side f.
1397          B^* is the adjoint operator of B is the given inner product.
1398    
1399          """
1400          def __init__(self,**kwargs):
1401            self.setTolerance()
1402            self.setToleranceReductionFactor()
1403    
1404          def initialize(self):
1405            """
1406            initialize the problem (overwrite)
1407            """
1408            pass
1409          def B(self,v):
1410             """
1411             returns Bv (overwrite)
1412             @rtype: equal to the type of p
1413    
1414             @note: boundary conditions on p should be zero!
1415             """
1416             pass
1417    
1418          def inner(self,p0,p1):
1419             """
1420             returns inner product of two element p0 and p1  (overwrite)
1421            
1422             @type p0: equal to the type of p
1423             @type p1: equal to the type of p
1424             @rtype: C{float}
1425    
1426             @rtype: equal to the type of p
1427             """
1428             pass
1429    
1430          def solve_A(self,u,p):
1431             """
1432             solves Av=f-Au-B^*p with accuracy self.getReducedTolerance() (overwrite)
1433    
1434             @rtype: equal to the type of v
1435             @note: boundary conditions on v should be zero!
1436             """
1437             pass
1438    
1439          def solve_prec(self,p):
1440             """
1441             provides a preconditioner for BA^{-1}B^* with accuracy self.getReducedTolerance() (overwrite)
1442    
1443             @rtype: equal to the type of p
1444             """
1445             pass
1446    
1447          def stoppingcriterium(self,Bv,v,p):
1448             """
1449             returns a True if iteration is terminated. (overwrite)
1450    
1451             @rtype: C{bool}
1452             """
1453             pass
1454                
1455          def __inner(self,p,r):
1456             return self.inner(p,r[1])
1457    
1458          def __inner_p(self,p1,p2):
1459             return self.inner(p1,p2)
1460          
1461          def __inner_a(self,a1,a2):
1462             return self.inner_a(a1,a2)
1463    
1464          def __inner_a1(self,a1,a2):
1465             return self.inner(a1[1],a2[1])
1466    
1467          def __stoppingcriterium(self,norm_r,r,p):
1468              return self.stoppingcriterium(r[1],r[0],p)
1469    
1470          def __stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None):
1471              return self.stoppingcriterium2(norm_r,norm_b,solver,TOL)
1472    
1473          def setTolerance(self,tolerance=1.e-8):
1474                  self.__tol=tolerance
1475          def getTolerance(self):
1476                  return self.__tol
1477          def setToleranceReductionFactor(self,reduction=0.01):
1478                  self.__reduction=reduction
1479          def getSubProblemTolerance(self):
1480                  return self.__reduction*self.getTolerance()
1481    
1482          def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='PCG',iter_restart=20):
1483                  """
1484                  solves the saddle point problem using initial guesses v and p.
1485    
1486                  @param max_iter: maximum number of iteration steps.
1487                  """
1488                  self.verbose=verbose
1489                  self.show_details=show_details and self.verbose
1490    
1491                  # assume p is known: then v=A^-1(f-B^*p)
1492                  # which leads to BA^-1B^*p = BA^-1f  
1493    
1494              # Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask)      
1495              self.__z=v+self.solve_A(v,p*0)
1496                  Bz=self.B(self.__z)
1497                  #
1498              #   solve BA^-1B^*p = Bz
1499                  #
1500                  #
1501                  #
1502                  self.iter=0
1503              if solver=='GMRES':      
1504                    if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter
1505                    p=GMRES(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.,iter_restart=iter_restart)
1506                    # solve Au=f-B^*p
1507                    #       A(u-v)=f-B^*p-Av
1508                    #       u=v+(u-v)
1509            u=v+self.solve_A(v,p)
1510    
1511              if solver=='TFQMR':      
1512                    if self.verbose: print "enter TFQMR method (iter_max=%s)"%max_iter
1513                    p=TFQMR(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)
1514                    # solve Au=f-B^*p
1515                    #       A(u-v)=f-B^*p-Av
1516                    #       u=v+(u-v)
1517            u=v+self.solve_A(v,p)
1518    
1519              if solver=='MINRES':      
1520                    if self.verbose: print "enter MINRES method (iter_max=%s)"%max_iter
1521                    p=MINRES(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)
1522                    # solve Au=f-B^*p
1523                    #       A(u-v)=f-B^*p-Av
1524                    #       u=v+(u-v)
1525            u=v+self.solve_A(v,p)
1526                  
1527              if solver=='GMRESC':      
1528                    if self.verbose: print "enter GMRES coupled method (iter_max=%s)"%max_iter
1529                    p0=self.solve_prec1(Bz)
1530                #alfa=(1/self.vol)*util.integrate(util.interpolate(p,escript.Function(self.domain)))
1531                    #p-=alfa
1532                    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)
1533                    #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())
1534    
1535                    # solve Au=f-B^*p
1536                    #       A(u-v)=f-B^*p-Av
1537                    #       u=v+(u-v)
1538                p=x[1]
1539            u=v+self.solve_A(v,p)      
1540            #p=x[1]
1541            #u=x[0]
1542    
1543                  if solver=='PCG':
1544                    #   note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv
1545                    #
1546                    #   with                    Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask)
1547                    #                           A(v-z)= f -Az - B^*p (v-z=0 on fixed_u_mask)
1548                    if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter
1549                    p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p)
1550                u=r[0]  
1551                    # print "DIFF=",util.integrate(p)
1552    
1553                  # print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)
1554    
1555              return u,p
1556    
1557          def __Msolve(self,r):
1558              return self.solve_prec(r[1])
1559    
1560          def __Msolve2(self,r):
1561              return self.solve_prec(r*1)
1562    
1563          def __Mempty(self,r):
1564              return r
1565    
1566    
1567          def __Aprod(self,p):
1568              # return BA^-1B*p
1569              #solve Av =B^*p as Av =f-Az-B^*(-p)
1570              v=self.solve_A(self.__z,-p)
1571              return ArithmeticTuple(v, self.B(v))
1572    
1573          def __Anext(self,x):
1574              # return next v,p
1575              #solve Av =-B^*p as Av =f-Az-B^*p
1576    
1577          pc=x[1]
1578              v=self.solve_A(self.__z,-pc)
1579          p=self.solve_prec1(self.B(v))
1580    
1581              return ArithmeticTuple(v,p)
1582    
1583    
1584          def __Aprod2(self,p):
1585              # return BA^-1B*p
1586              #solve Av =B^*p as Av =f-Az-B^*(-p)
1587          v=self.solve_A(self.__z,-p)
1588              return self.B(v)
1589    
1590          def __Aprod_Newton(self,p):
1591              # return BA^-1B*p - Bz
1592              #solve Av =-B^*p as Av =f-Az-B^*p
1593          v=self.solve_A(self.__z,-p)
1594              return self.B(v-self.__z)
1595    
1596          def __Aprod_Newton2(self,x):
1597              # return BA^-1B*p - Bz
1598              #solve Av =-B^*p as Av =f-Az-B^*p
1599              pc=x[1]
1600          v=self.solve_A(self.__z,-pc)
1601              p=self.solve_prec1(self.B(v-self.__z))
1602              return ArithmeticTuple(v,p)
1603    
1604    
1605    def MaskFromBoundaryTag(domain,*tags):
1606       """
1607       creates a mask on the Solution(domain) function space which one for samples
1608       that touch the boundary tagged by tags.
1609    
1610       usage: m=MaskFromBoundaryTag(domain,"left", "right")
1611    
1612       @param domain: a given domain
1613       @type domain: L{escript.Domain}
1614       @param tags: boundray tags
1615       @type tags: C{str}
1616       @return: a mask which marks samples that are touching the boundary tagged by any of the given tags.
1617       @rtype: L{escript.Data} of rank 0
1618       """
1619       pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1620       d=escript.Scalar(0.,escript.FunctionOnBoundary(domain))
1621       for t in tags: d.setTaggedValue(t,1.)
1622       pde.setValue(y=d)
1623       return util.whereNonZero(pde.getRightHandSide())
1624    #============================================================================================================================================
1625    class SaddlePointProblem(object):
1626       """
1627       This implements a solver for a saddlepoint problem
1628    
1629       M{f(u,p)=0}
1630       M{g(u)=0}
1631    
1632       for u and p. The problem is solved with an inexact Uszawa scheme for p:
1633    
1634       M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})}
1635       M{Q_g (p^{k+1}-p^{k}) =   g(u^{k+1})}
1636    
1637       where Q_f is an approximation of the Jacobiean A_f of f with respect to u  and Q_f is an approximation of
1638       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'
1639       Q_g can be difficult, non-linear conjugate gradient method is applied to solve for p, so Q_g plays
1640       in fact the role of a preconditioner.
1641       """
1642       def __init__(self,verbose=False,*args):
1643           """
1644           initializes the problem
1645    
1646           @param verbose: switches on the printing out some information
1647           @type verbose: C{bool}
1648           @note: this method may be overwritten by a particular saddle point problem
1649           """
1650           print "SaddlePointProblem should not be used anymore!"
1651           if not isinstance(verbose,bool):
1652                raise TypeError("verbose needs to be of type bool.")
1653           self.__verbose=verbose
1654           self.relaxation=1.
1655           DeprecationWarning("SaddlePointProblem should not be used anymore.")
1656    
1657       def trace(self,text):
1658           """
1659           prints text if verbose has been set
1660    
1661           @param text: a text message
1662           @type text: C{str}
1663           """
1664           if self.__verbose: print "%s: %s"%(str(self),text)
1665    
1666       def solve_f(self,u,p,tol=1.e-8):
1667           """
1668           solves
1669    
1670           A_f du = f(u,p)
1671    
1672           with tolerance C{tol} and return du. A_f is Jacobiean of f with respect to u.
1673    
1674           @param u: current approximation of u
1675           @type u: L{escript.Data}
1676           @param p: current approximation of p
1677           @type p: L{escript.Data}
1678           @param tol: tolerance expected for du
1679           @type tol: C{float}
1680           @return: increment du
1681           @rtype: L{escript.Data}
1682           @note: this method has to be overwritten by a particular saddle point problem
1683           """
1684           pass
1685    
1686       def solve_g(self,u,tol=1.e-8):
1687           """
1688           solves
1689    
1690           Q_g dp = g(u)
1691    
1692           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.
1693    
1694           @param u: current approximation of u
1695           @type u: L{escript.Data}
1696           @param tol: tolerance expected for dp
1697           @type tol: C{float}
1698           @return: increment dp
1699           @rtype: L{escript.Data}
1700           @note: this method has to be overwritten by a particular saddle point problem
1701           """
1702           pass
1703    
1704       def inner(self,p0,p1):
1705           """
1706           inner product of p0 and p1 approximating p. Typically this returns integrate(p0*p1)
1707           @return: inner product of p0 and p1
1708           @rtype: C{float}
1709           """
1710           pass
1711    
1712       subiter_max=3
1713       def solve(self,u0,p0,tolerance=1.e-6,tolerance_u=None,iter_max=100,accepted_reduction=0.995,relaxation=None):
1714            """
1715            runs the solver
1716    
1717            @param u0: initial guess for C{u}
1718            @type u0: L{esys.escript.Data}
1719            @param p0: initial guess for C{p}
1720            @type p0: L{esys.escript.Data}
1721            @param tolerance: tolerance for relative error in C{u} and C{p}
1722            @type tolerance: positive C{float}
1723            @param tolerance_u: tolerance for relative error in C{u} if different from C{tolerance}
1724            @type tolerance_u: positive C{float}
1725            @param iter_max: maximum number of iteration steps.
1726            @type iter_max: C{int}
1727            @param accepted_reduction: if the norm  g cannot be reduced by C{accepted_reduction} backtracking to adapt the
1728                                       relaxation factor. If C{accepted_reduction=None} no backtracking is used.
1729            @type accepted_reduction: positive C{float} or C{None}
1730            @param relaxation: initial relaxation factor. If C{relaxation==None}, the last relaxation factor is used.
1731            @type relaxation: C{float} or C{None}
1732            """
1733            tol=1.e-2
1734            if tolerance_u==None: tolerance_u=tolerance
1735            if not relaxation==None: self.relaxation=relaxation
1736            if accepted_reduction ==None:
1737                  angle_limit=0.
1738            elif accepted_reduction>=1.:
1739                  angle_limit=0.
1740            else:
1741                  angle_limit=util.sqrt(1-accepted_reduction**2)
1742            self.iter=0
1743            u=u0
1744            p=p0
1745            #
1746            #   initialize things:
1747            #
1748            converged=False
1749            #
1750            #  start loop:
1751            #
1752            #  initial search direction is g
1753            #
1754            while not converged :
1755                if self.iter>iter_max:
1756                    raise ArithmeticError("no convergence after %s steps."%self.iter)
1757                f_new=self.solve_f(u,p,tol)
1758                norm_f_new = util.Lsup(f_new)
1759                u_new=u-f_new
1760                g_new=self.solve_g(u_new,tol)
1761                self.iter+=1
1762                norm_g_new = util.sqrt(self.inner(g_new,g_new))
1763                if norm_f_new==0. and norm_g_new==0.: return u, p
1764                if self.iter>1 and not accepted_reduction==None:
1765                   #
1766                   #   did we manage to reduce the norm of G? I
1767                   #   if not we start a backtracking procedure
1768                   #
1769                   # print "new/old norm = ",norm_g_new, norm_g, norm_g_new/norm_g
1770                   if norm_g_new > accepted_reduction * norm_g:
1771                      sub_iter=0
1772                      s=self.relaxation
1773                      d=g
1774                      g_last=g
1775                      self.trace("    start substepping: f = %s, g = %s, relaxation = %s."%(norm_f_new, norm_g_new, s))
1776                      while sub_iter < self.subiter_max and  norm_g_new > accepted_reduction * norm_g:
1777                         dg= g_new-g_last
1778                         norm_dg=abs(util.sqrt(self.inner(dg,dg))/self.relaxation)
1779                         rad=self.inner(g_new,dg)/self.relaxation
1780                         # print "   ",sub_iter,": rad, norm_dg:",abs(rad), norm_dg*norm_g_new * angle_limit
1781                         # print "   ",sub_iter,": rad, norm_dg:",rad, norm_dg, norm_g_new, norm_g
1782                         if abs(rad) < norm_dg*norm_g_new * angle_limit:
1783                             if sub_iter>0: self.trace("    no further improvements expected from backtracking.")
1784                             break
1785                         r=self.relaxation
1786                         self.relaxation= - rad/norm_dg**2
1787                         s+=self.relaxation
1788                         #####
1789                         # a=g_new+self.relaxation*dg/r
1790                         # print "predicted new norm = ",util.sqrt(self.inner(a,a)),util.sqrt(self.inner(g_new,g_new)), self.relaxation
1791                         #####
1792                         g_last=g_new
1793                         p+=self.relaxation*d
1794                         f_new=self.solve_f(u,p,tol)
1795                         u_new=u-f_new
1796                         g_new=self.solve_g(u_new,tol)
1797                         self.iter+=1
1798                         norm_f_new = util.Lsup(f_new)
1799                         norm_g_new = util.sqrt(self.inner(g_new,g_new))
1800                         # print "   ",sub_iter," new g norm",norm_g_new
1801                         self.trace("    %s th sub-step: f = %s, g = %s, relaxation = %s."%(sub_iter, norm_f_new, norm_g_new, s))
1802                         #
1803                         #   can we expect reduction of g?
1804                         #
1805                         # u_last=u_new
1806                         sub_iter+=1
1807                      self.relaxation=s
1808                #
1809                #  check for convergence:
1810                #
1811                norm_u_new = util.Lsup(u_new)
1812                p_new=p+self.relaxation*g_new
1813                norm_p_new = util.sqrt(self.inner(p_new,p_new))
1814                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))
1815    
1816                if self.iter>1:
1817                   dg2=g_new-g
1818                   df2=f_new-f
1819                   norm_dg2=util.sqrt(self.inner(dg2,dg2))
1820                   norm_df2=util.Lsup(df2)
1821                   # print norm_g_new, norm_g, norm_dg, norm_p, tolerance
1822                   tol_eq_g=tolerance*norm_dg2/(norm_g*abs(self.relaxation))*norm_p_new
1823                   tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new
1824                   if norm_g_new <= tol_eq_g and norm_f_new <= tol_eq_f:
1825                       converged=True
1826                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
1827            self.trace("convergence after %s steps."%self.iter)
1828            return u,p

Legend:
Removed from v.155  
changed lines
  Added in v.1956

  ViewVC Help
Powered by ViewVC 1.1.26