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

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

  ViewVC Help
Powered by ViewVC 1.1.26