/[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 396 by gross, Wed Dec 21 05:08:25 2005 UTC revision 1552 by gross, Thu May 8 08:52:41 2008 UTC
# Line 1  Line 1 
1    #
2  # $Id$  # $Id$
3    #
4    #######################################################
5    #
6    #           Copyright 2003-2007 by ACceSS MNRF
7    #       Copyright 2007 by University of Queensland
8    #
9    #                http://esscc.uq.edu.au
10    #        Primary Business: Queensland, Australia
11    #  Licensed under the Open Software License version 3.0
12    #     http://www.opensource.org/licenses/osl-3.0.php
13    #
14    #######################################################
15    #
16    
17  """  """
18  Provides some tools related to PDEs.  Provides some tools related to PDEs.
# Line 7  Currently includes: Line 21  Currently includes:
21      - Projector - to project a discontinuous      - Projector - to project a discontinuous
22      - Locator - to trace values in data objects at a certain location      - Locator - to trace values in data objects at a certain location
23      - TimeIntegrationManager - to handel extraplotion in time      - TimeIntegrationManager - to handel extraplotion in time
24        - SaddlePointProblem - solver for Saddle point problems using the inexact uszawa scheme
25    
26    @var __author__: name of author
27    @var __copyright__: copyrights
28    @var __license__: licence agreement
29    @var __url__: url entry point on documentation
30    @var __version__: version
31    @var __date__: date of the version
32  """  """
33    
34    __author__="Lutz Gross, l.gross@uq.edu.au"
35    __copyright__="""  Copyright (c) 2006 by ACcESS MNRF
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:  class TimeIntegrationManager:
57    """    """
58    a simple mechanism to manage time dependend values.    a simple mechanism to manage time dependend values.
59    
60    typical usage is:    typical usage is::
61    
62    dt=0.1 # time increment       dt=0.1 # time increment
63    tm=TimeIntegrationManager(inital_value,p=1)       tm=TimeIntegrationManager(inital_value,p=1)
64    while t<1.       while t<1.
65        v_guess=tm.extrapolate(dt) # extrapolate to t+dt           v_guess=tm.extrapolate(dt) # extrapolate to t+dt
66        v=...           v=...
67        tm.checkin(dt,v)           tm.checkin(dt,v)
68        t+=dt           t+=dt
69    
70    @remark: currently only p=1 is supported.    @note: currently only p=1 is supported.
71    """    """
72    def __init__(self,*inital_values,**kwargs):    def __init__(self,*inital_values,**kwargs):
73       """       """
# Line 50  class TimeIntegrationManager: Line 89  class TimeIntegrationManager:
89    def getTime(self):    def getTime(self):
90        return self.__t        return self.__t
91    def getValue(self):    def getValue(self):
92        return self.__v_mem[0]        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):    def checkin(self,dt,*values):
99        """        """
100        adds new values to the manager. the p+1 last value get lost        adds new values to the manager. the p+1 last value get lost
# Line 108  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 145  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 161  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 180  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 206  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       while not (stoppingcriterium(rho,norm_b) or iter==iter_restart-1):
622    
623        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
624    
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    
641    # Reorthogonalize if needed
642        if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
643         for j in range(iter+1):
644            hr=bilinearform(v[j],v[iter+1])
645                h[j][iter]=h[j][iter]+hr #vhat
646                v[iter+1] -= hr*v[j]
647    
648         v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
649         h[iter+1][iter]=v_norm2
650    
651    #   watch out for happy breakdown
652            if not v_norm2 == 0:
653             v[iter+1]=v[iter+1]/h[iter+1][iter]
654    
655    #   Form and store the information for the new Givens rotation
656        if iter > 0 :
657            hhat=numarray.zeros(iter+1,numarray.Float64)
658            for i in range(iter+1) : hhat[i]=h[i][iter]
659            hhat=givapp(c[0:iter],s[0:iter],hhat);
660                for i in range(iter+1) : h[i][iter]=hhat[i]
661    
662        mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])
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            rho=abs(g[iter+1])
672        iter+=1
673    
674    # At this point either iter > iter_max or rho < tol.
675    # It's time to compute x and leave.        
676    
677       if iter > 0 :
678         y=numarray.zeros(iter,numarray.Float64)    
679         y[iter-1] = g[iter-1] / h[iter-1][iter-1]
680         if iter > 1 :  
681            i=iter-2  
682            while i>=0 :
683              y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]
684              i=i-1
685         xhat=v[iter-1]*y[iter-1]
686         for i in range(iter-1):
687        xhat += v[i]*y[i]
688       else : xhat=v[0]
689        
690       x += xhat
691       if iter<iter_restart-1:
692          stopped=True
693       else:
694          stopped=False
695    
696       return x,stopped
697    
698    
699    ######################################################
700    def dirder(x, w, bilinearform, Aprod, Msolve, f0, b ):
701    ######################################################
702    
703    # DIRDER estimates the directional derivative of a function F.
704    
705    
706    # Hardwired difference increment.
707    #
708      epsnew = 1.0e-07
709    #
710    #  Scale the step.
711    #
712      norm_w=math.sqrt(bilinearform(w,w))
713      if norm_w== 0.0:
714        return x/x
715    
716      epsnew = epsnew / norm_w
717    
718      if norm_w > 0.0:
719        epsnew = epsnew * math.sqrt(bilinearform(x,x))
720    #
721    #  DEL and F1 could share the same space if storage
722    #  is more important than clarity.
723    #
724    
725      DEL = x + epsnew * w
726      f1 = -Msolve(Aprod(DEL))
727      z = ( f1 - f0 ) / epsnew
728      return z
729    
730    ######################################################
731    def FDGMRES(f0, Aprod, Msolve, bilinearform, stoppingcriterium, xc=None, x=None, iter_max=100, iter_restart=20,TOL=None):
732    ######################################################
733       b=-f0
734       b_dot_b = bilinearform(b, b)
735       if b_dot_b<0: raise NegativeNorm,"negative norm."
736       norm_b=math.sqrt(b_dot_b)
737    
738       r=b
739    
740       if x==None:
741          x=0*f0
742       else:
743          r=-dirder(xc,x,bilinearform,Aprod,Msolve,f0,b)-f0  
744          
745       r_dot_r = bilinearform(r, r)
746       if r_dot_r<0: raise NegativeNorm,"negative norm."
747      
748       h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)
749       c=numarray.zeros(iter_restart,numarray.Float64)
750       s=numarray.zeros(iter_restart,numarray.Float64)
751       g=numarray.zeros(iter_restart,numarray.Float64)
752       v=[]
753    
754       rho=math.sqrt(r_dot_r)
755      
756       v.append(r/rho)
757       g[0]=rho
758       iter=0
759    
760       while not (stoppingcriterium(rho,norm_b,solver="FDGMRES",TOL=TOL) or iter==iter_restart-1):
761    
762        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
763    
764        
765            p=dirder(xc, v[iter], bilinearform,Aprod,Msolve,f0,b)
766    
767        v.append(p)
768    
769        v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
770    
771    # Modified Gram-Schmidt
772        for j in range(iter+1):
773          h[j][iter]=bilinearform(v[j],v[iter+1])  
774          v[iter+1]+=(-1.)*h[j][iter]*v[j]
775          
776        h[iter+1][iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
777        v_norm2=h[iter+1][iter]
778    
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=bilinearform(v[j],v[iter+1])
784                h[j][iter]=h[j][iter]+hr #vhat
785                v[iter+1] +=(-1.)*hr*v[j]
786    
787         v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
788         h[iter+1][iter]=v_norm2
789    
790    #   watch out for happy breakdown
791            if v_norm2 != 0:
792             v[iter+1]=v[iter+1]/h[iter+1][iter]
793    
794    #   Form and store the information for the new Givens rotation
795        if iter > 0 :
796            hhat=[]
797            for i in range(iter+1) : hhat.append(h[i][iter])
798            hhat=givapp(c[0:iter],s[0:iter],hhat);
799                for i in range(iter+1) : h[i][iter]=hhat[i]
800    
801        mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])
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       if iter > 0 :
814         y=numarray.zeros(iter,numarray.Float64)    
815         y[iter-1] = g[iter-1] / h[iter-1][iter-1]
816         if iter > 1 :  
817            i=iter-2  
818            while i>=0 :
819              y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]
820              i=i-1
821         xhat=v[iter-1]*y[iter-1]
822         for i in range(iter-1):
823        xhat += v[i]*y[i]
824       else : xhat=v[0]
825        
826       x += xhat
827       if iter<iter_restart-1:
828          stopped=True
829       else:
830          stopped=False
831    
832       return x,stopped
833    
834    #################################################
835    def MINRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
836    #################################################
837        #
838        #  minres solves the system of linear equations Ax = b
839        #  where A is a symmetric matrix (possibly indefinite or singular)
840        #  and b is a given vector.
841        #  
842        #  "A" may be a dense or sparse matrix (preferably sparse!)
843        #  or the name of a function such that
844        #               y = A(x)
845        #  returns the product y = Ax for any given vector x.
846        #
847        #  "M" defines a positive-definite preconditioner M = C C'.
848        #  "M" may be a dense or sparse matrix (preferably sparse!)
849        #  or the name of a function such that
850        #  solves the system My = x for any given vector x.
851        #
852        #
853        
854        #------------------------------------------------------------------
855        # Set up y and v for the first Lanczos vector v1.
856        # y  =  beta1 P' v1,  where  P = C**(-1).
857        # v is really P' v1.
858        #------------------------------------------------------------------
859        if x==None:
860          x=0*b
861        else:
862          b += (-1)*Aprod(x)
863    
864        r1    = b
865        y = Msolve(b)
866        beta1 = bilinearform(y,b)
867    
868        if beta1< 0: raise NegativeNorm,"negative norm."
869    
870        #  If b = 0 exactly, stop with x = 0.
871        if beta1==0: return x*0.
872    
873        if beta1> 0:
874          beta1  = math.sqrt(beta1)      
875    
876        #------------------------------------------------------------------
877        # Initialize quantities.
878        # ------------------------------------------------------------------
879        iter   = 0
880        Anorm = 0
881        ynorm = 0
882        oldb   = 0
883        beta   = beta1
884        dbar   = 0
885        epsln  = 0
886        phibar = beta1
887        rhs1   = beta1
888        rhs2   = 0
889        rnorm  = phibar
890        tnorm2 = 0
891        ynorm2 = 0
892        cs     = -1
893        sn     = 0
894        w      = b*0.
895        w2     = b*0.
896        r2     = r1
897        eps    = 0.0001
898    
899        #---------------------------------------------------------------------
900        # Main iteration loop.
901        # --------------------------------------------------------------------
902        while not stoppingcriterium(rnorm,Anorm*ynorm,'MINRES'):    #  checks ||r|| < (||A|| ||x||) * TOL
903    
904        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
905            iter    = iter  +  1
906    
907            #-----------------------------------------------------------------
908            # Obtain quantities for the next Lanczos vector vk+1, k = 1, 2,...
909            # The general iteration is similar to the case k = 1 with v0 = 0:
910            #
911            #   p1      = Operator * v1  -  beta1 * v0,
912            #   alpha1  = v1'p1,
913            #   q2      = p2  -  alpha1 * v1,
914            #   beta2^2 = q2'q2,
915            #   v2      = (1/beta2) q2.
916            #
917            # Again, y = betak P vk,  where  P = C**(-1).
918            #-----------------------------------------------------------------
919            s = 1/beta                 # Normalize previous vector (in y).
920            v = s*y                    # v = vk if P = I
921        
922            y      = Aprod(v)
923        
924            if iter >= 2:
925              y = y - (beta/oldb)*r1
926    
927            alfa   = bilinearform(v,y)              # alphak
928            y      += (- alfa/beta)*r2
929            r1     = r2
930            r2     = y
931            y = Msolve(r2)
932            oldb   = beta                           # oldb = betak
933            beta   = bilinearform(y,r2)             # beta = betak+1^2
934            if beta < 0: raise NegativeNorm,"negative norm."
935    
936            beta   = math.sqrt( beta )
937            tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta
938            
939            if iter==1:                 # Initialize a few things.
940              gmax   = abs( alfa )      # alpha1
941              gmin   = gmax             # alpha1
942    
943            # Apply previous rotation Qk-1 to get
944            #   [deltak epslnk+1] = [cs  sn][dbark    0   ]
945            #   [gbar k dbar k+1]   [sn -cs][alfak betak+1].
946        
947            oldeps = epsln
948            delta  = cs * dbar  +  sn * alfa  # delta1 = 0         deltak
949            gbar   = sn * dbar  -  cs * alfa  # gbar 1 = alfa1     gbar k
950            epsln  =               sn * beta  # epsln2 = 0         epslnk+1
951            dbar   =            -  cs * beta  # dbar 2 = beta2     dbar k+1
952    
953            # Compute the next plane rotation Qk
954    
955            gamma  = math.sqrt(gbar*gbar+beta*beta)  # gammak
956            gamma  = max(gamma,eps)
957            cs     = gbar / gamma             # ck
958            sn     = beta / gamma             # sk
959            phi    = cs * phibar              # phik
960            phibar = sn * phibar              # phibark+1
961    
962            # Update  x.
963    
964            denom = 1/gamma
965            w1    = w2
966            w2    = w
967            w     = (v - oldeps*w1 - delta*w2) * denom
968            x     +=  phi*w
969    
970            # Go round again.
971    
972            gmax   = max(gmax,gamma)
973            gmin   = min(gmin,gamma)
974            z      = rhs1 / gamma
975            ynorm2 = z*z  +  ynorm2
976            rhs1   = rhs2 -  delta*z
977            rhs2   =      -  epsln*z
978    
979            # Estimate various norms and test for convergence.
980    
981            Anorm  = math.sqrt( tnorm2 )
982            ynorm  = math.sqrt( ynorm2 )
983    
984            rnorm  = phibar
985    
986        return x
987    
988    ######################################    
989    def NewtonGMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium,x=None, iter_max=100,iter_restart=20,atol=1.e-2,rtol=1.e-4):
990    #####################################
991        gamma=.9
992        lmaxit=100
993        etamax=.5
994    
995        n = 1 #len(x)
996        iter=0
997        
998        # evaluate f at the initial iterate
999        # compute the stop tolerance
1000        #
1001        r=b
1002        if x==None:
1003          x=0*b
1004        else:
1005          r += (-1)*Aprod(x)
1006    
1007        f0=-Msolve(r)
1008        fnrm=math.sqrt(bilinearform(f0,f0))/math.sqrt(n)
1009        fnrmo=1
1010        stop_tol=atol + rtol*fnrm
1011        #
1012        # main iteration loop
1013        #
1014        while not stoppingcriterium(fnrm*1,stop_tol,'NewtonGMRES',TOL=1.):
1015    
1016                if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
1017            #
1018            # keep track of the ratio (rat = fnrm/frnmo)
1019            # of successive residual norms and
1020            # the iteration counter (iter)
1021            #
1022            #rat=fnrm/fnrmo
1023            fnrmo=fnrm
1024            iter+=1
1025            #
1026                # compute the step using a GMRES(m) routine especially designed for this purpose
1027            #
1028                initer=0
1029                while True:
1030                   xc,stopped=FDGMRES(f0*1, Aprod, Msolve, bilinearform, stoppingcriterium, xc=x*1, iter_max=lmaxit-initer, iter_restart=iter_restart, TOL=etamax)
1031                   if stopped: break
1032                   initer+=iter_restart
1033            x+=xc
1034            f0=-Msolve(Aprod(x))
1035            fnrm=math.sqrt(bilinearform(f0,f0))/math.sqrt(n)
1036            rat=fnrm/fnrmo
1037    
1038    
1039        #   adjust eta
1040        #
1041            if etamax > 0:
1042                etaold=etamax
1043                etanew=gamma*rat*rat
1044                if gamma*etaold*etaold > .1 :
1045                    etanew=max(etanew,gamma*etaold*etaold)
1046                
1047                etamax=min(etanew,etamax)
1048                etamax=max(etamax,.5*stop_tol/fnrm)
1049    
1050        return x
1051    
1052    def TFQMR(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
1053    
1054    # TFQMR solver for linear systems
1055    #
1056    #
1057    # initialization
1058    #
1059      errtol = math.sqrt(bilinearform(b,b))
1060      norm_b=errtol
1061      kmax  = iter_max
1062      error = []
1063    
1064      if math.sqrt(bilinearform(x,x)) != 0.0:
1065        r = b - Aprod(x)
1066      else:
1067        r = b
1068    
1069      r=Msolve(r)
1070    
1071      u1=0
1072      u2=0
1073      y1=0
1074      y2=0
1075    
1076      w = r
1077      y1 = r
1078      iter = 0
1079      d = 0
1080      
1081      v = Msolve(Aprod(y1))
1082      u1 = v
1083      
1084      theta = 0.0;
1085      eta = 0.0;
1086      tau = math.sqrt(bilinearform(r,r))
1087      error = [ error, tau ]
1088      rho = tau * tau
1089      m=1
1090    #
1091    #  TFQMR iteration
1092    #
1093    #  while ( iter < kmax-1 ):
1094      
1095      while not stoppingcriterium(tau*math.sqrt ( m + 1 ),norm_b,'TFQMR'):
1096        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
1097    
1098        sigma = bilinearform(r,v)
1099    
1100        if ( sigma == 0.0 ):
1101          raise 'TFQMR breakdown, sigma=0'
1102        
1103    
1104        alpha = rho / sigma
1105    
1106        for j in range(2):
1107    #
1108    #   Compute y2 and u2 only if you have to
1109    #
1110          if ( j == 1 ):
1111            y2 = y1 - alpha * v
1112            u2 = Msolve(Aprod(y2))
1113    
1114          m = 2 * (iter+1) - 2 + (j+1)
1115          if j==0:
1116             w = w - alpha * u1
1117             d = y1 + ( theta * theta * eta / alpha ) * d
1118          if j==1:
1119             w = w - alpha * u2
1120             d = y2 + ( theta * theta * eta / alpha ) * d
1121    
1122          theta = math.sqrt(bilinearform(w,w))/ tau
1123          c = 1.0 / math.sqrt ( 1.0 + theta * theta )
1124          tau = tau * theta * c
1125          eta = c * c * alpha
1126          x = x + eta * d
1127    #
1128    #  Try to terminate the iteration at each pass through the loop
1129    #
1130         # if ( tau * math.sqrt ( m + 1 ) <= errtol ):
1131         #   error = [ error, tau ]
1132         #   total_iters = iter
1133         #   break
1134          
1135    
1136        if ( rho == 0.0 ):
1137          raise 'TFQMR breakdown, rho=0'
1138        
1139    
1140        rhon = bilinearform(r,w)
1141        beta = rhon / rho;
1142        rho = rhon;
1143        y1 = w + beta * y2;
1144        u1 = Msolve(Aprod(y1))
1145        v = u1 + beta * ( u2 + beta * v )
1146        error = [ error, tau ]
1147        total_iters = iter
1148        
1149        iter = iter + 1
1150    
1151      return x
1152    
1153    
1154    #############################################
1155    
1156    class ArithmeticTuple(object):
1157       """
1158       tuple supporting inplace update x+=y and scaling x=a*y where x,y is an ArithmeticTuple and a is a float.
1159    
1160       example of usage:
1161    
1162       from esys.escript import Data
1163       from numarray import array
1164       a=Data(...)
1165       b=array([1.,4.])
1166       x=ArithmeticTuple(a,b)
1167       y=5.*x
1168    
1169       """
1170       def __init__(self,*args):
1171           """
1172           initialize object with elements args.
1173    
1174           @param args: tuple of object that support implace add (x+=y) and scaling (x=a*y)
1175           """
1176           self.__items=list(args)
1177    
1178       def __len__(self):
1179           """
1180           number of items
1181    
1182           @return: number of items
1183           @rtype: C{int}
1184           """
1185           return len(self.__items)
1186    
1187       def __getitem__(self,index):
1188           """
1189           get an item
1190    
1191           @param index: item to be returned
1192           @type index: C{int}
1193           @return: item with index C{index}
1194           """
1195           return self.__items.__getitem__(index)
1196    
1197       def __mul__(self,other):
1198           """
1199           scaling from the right
1200    
1201           @param other: scaling factor
1202           @type other: C{float}
1203           @return: itemwise self*other
1204           @rtype: L{ArithmeticTuple}
1205           """
1206           out=[]
1207           for i in range(len(self)):
1208               out.append(self[i]*other)
1209           return ArithmeticTuple(*tuple(out))
1210    
1211       def __rmul__(self,other):
1212           """
1213           scaling from the left
1214    
1215           @param other: scaling factor
1216           @type other: C{float}
1217           @return: itemwise other*self
1218           @rtype: L{ArithmeticTuple}
1219           """
1220           out=[]
1221           for i in range(len(self)):
1222               out.append(other*self[i])
1223           return ArithmeticTuple(*tuple(out))
1224    
1225    #########################
1226    # Added by Artak
1227    #########################
1228       def __div__(self,other):
1229           """
1230           dividing from the right
1231    
1232           @param other: scaling factor
1233           @type other: C{float}
1234           @return: itemwise self/other
1235           @rtype: L{ArithmeticTuple}
1236           """
1237           out=[]
1238           for i in range(len(self)):
1239               out.append(self[i]*(1./other))
1240           return ArithmeticTuple(*tuple(out))
1241    
1242       def __rdiv__(self,other):
1243           """
1244           dividing from the left
1245    
1246           @param other: scaling factor
1247           @type other: C{float}
1248           @return: itemwise other/self
1249           @rtype: L{ArithmeticTuple}
1250           """
1251           out=[]
1252           for i in range(len(self)):
1253               out.append(other/self[i])
1254           return ArithmeticTuple(*tuple(out))
1255      
1256    ##########################################33
1257    
1258       def __iadd__(self,other):
1259           """
1260           in-place add of other to self
1261    
1262           @param other: increment
1263           @type other: C{ArithmeticTuple}
1264           """
1265           if len(self) != len(other):
1266               raise ValueError,"tuple length must match."
1267           for i in range(len(self)):
1268               self.__items[i]+=other[i]
1269           return self
1270    
1271       def __add__(self,other):
1272           """
1273           add other to self
1274    
1275           @param other: increment
1276           @type other: C{ArithmeticTuple}
1277           """
1278           if not isinstance(other,float):
1279            if len(self) != len(other):
1280               raise ValueError,"tuple length must match."
1281            for i in range(len(self)):
1282               self.__items[i]+=other[i]
1283           else:
1284            for i in range(len(self)):
1285               self.__items[i]+=other
1286    
1287           return self
1288    
1289       def __sub__(self,other):
1290           """
1291           subtract other from self
1292    
1293           @param other: increment
1294           @type other: C{ArithmeticTuple}
1295           """
1296           if len(self) != len(other):
1297               raise ValueError,"tuple length must match."
1298           for i in range(len(self)):
1299               self.__items[i]-=other[i]
1300           return self
1301    
1302       def __neg__(self):
1303           """
1304           negate
1305    
1306           """
1307           for i in range(len(self)):
1308               self.__items[i]=-self.__items[i]
1309           return self
1310    
1311    
1312    class HomogeneousSaddlePointProblem(object):
1313          """
1314          This provides a framwork for solving homogeneous saddle point problem of the form
1315    
1316                 Av+B^*p=f
1317                 Bv    =0
1318    
1319          for the unknowns v and p and given operators A and B and given right hand side f.
1320          B^* is the adjoint operator of B is the given inner product.
1321    
1322          """
1323          def __init__(self,**kwargs):
1324            self.setTolerance()
1325            self.setToleranceReductionFactor()
1326    
1327          def initialize(self):
1328            """
1329            initialize the problem (overwrite)
1330            """
1331            pass
1332          def B(self,v):
1333             """
1334             returns Bv (overwrite)
1335             @rtype: equal to the type of p
1336    
1337             @note: boundary conditions on p should be zero!
1338             """
1339             pass
1340    
1341          def inner(self,p0,p1):
1342             """
1343             returns inner product of two element p0 and p1  (overwrite)
1344            
1345             @type p0: equal to the type of p
1346             @type p1: equal to the type of p
1347             @rtype: C{float}
1348    
1349             @rtype: equal to the type of p
1350             """
1351             pass
1352    
1353          def solve_A(self,u,p):
1354             """
1355             solves Av=f-Au-B^*p with accuracy self.getReducedTolerance() (overwrite)
1356    
1357             @rtype: equal to the type of v
1358             @note: boundary conditions on v should be zero!
1359             """
1360             pass
1361    
1362          def solve_prec(self,p):
1363             """
1364             provides a preconditioner for BA^{-1}B^* with accuracy self.getReducedTolerance() (overwrite)
1365    
1366             @rtype: equal to the type of p
1367             """
1368             pass
1369    
1370          def stoppingcriterium(self,Bv,v,p):
1371             """
1372             returns a True if iteration is terminated. (overwrite)
1373    
1374             @rtype: C{bool}
1375             """
1376             pass
1377                
1378          def __inner(self,p,r):
1379             return self.inner(p,r[1])
1380    
1381          def __inner_p(self,p1,p2):
1382             return self.inner(p1,p2)
1383          
1384          def __inner_a(self,a1,a2):
1385             return self.inner_a(a1,a2)
1386    
1387          def __inner_a1(self,a1,a2):
1388             return self.inner(a1[1],a2[1])
1389    
1390          def __stoppingcriterium(self,norm_r,r,p):
1391              return self.stoppingcriterium(r[1],r[0],p)
1392    
1393          def __stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None):
1394              return self.stoppingcriterium2(norm_r,norm_b,solver,TOL)
1395    
1396          def setTolerance(self,tolerance=1.e-8):
1397                  self.__tol=tolerance
1398          def getTolerance(self):
1399                  return self.__tol
1400          def setToleranceReductionFactor(self,reduction=0.01):
1401                  self.__reduction=reduction
1402          def getSubProblemTolerance(self):
1403                  return self.__reduction*self.getTolerance()
1404    
1405          def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='PCG',iter_restart=20):
1406                  """
1407                  solves the saddle point problem using initial guesses v and p.
1408    
1409                  @param max_iter: maximum number of iteration steps.
1410                  """
1411                  self.verbose=verbose
1412                  self.show_details=show_details and self.verbose
1413    
1414                  # assume p is known: then v=A^-1(f-B^*p)
1415                  # which leads to BA^-1B^*p = BA^-1f  
1416    
1417              # Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask)      
1418              self.__z=v+self.solve_A(v,p*0)
1419                  Bz=self.B(self.__z)
1420                  #
1421              #   solve BA^-1B^*p = Bz
1422                  #
1423                  #
1424                  #
1425                  self.iter=0
1426              if solver=='GMRES':      
1427                    if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter
1428                    p=GMRES(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.,iter_restart=iter_restart)
1429                    # solve Au=f-B^*p
1430                    #       A(u-v)=f-B^*p-Av
1431                    #       u=v+(u-v)
1432            u=v+self.solve_A(v,p)
1433    
1434              if solver=='NewtonGMRES':    
1435                    if self.verbose: print "enter NewtonGMRES method (iter_max=%s)"%max_iter
1436                    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())
1437                    # solve Au=f-B^*p
1438                    #       A(u-v)=f-B^*p-Av
1439                    #       u=v+(u-v)
1440            u=v+self.solve_A(v,p)
1441                    
1442    
1443              if solver=='TFQMR':      
1444                    if self.verbose: print "enter TFQMR method (iter_max=%s)"%max_iter
1445                    p=TFQMR(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)
1446                    # solve Au=f-B^*p
1447                    #       A(u-v)=f-B^*p-Av
1448                    #       u=v+(u-v)
1449            u=v+self.solve_A(v,p)
1450    
1451              if solver=='MINRES':      
1452                    if self.verbose: print "enter MINRES method (iter_max=%s)"%max_iter
1453                    p=MINRES(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)
1454                    # solve Au=f-B^*p
1455                    #       A(u-v)=f-B^*p-Av
1456                    #       u=v+(u-v)
1457            u=v+self.solve_A(v,p)
1458                  
1459              if solver=='GMRESC':      
1460                    if self.verbose: print "enter GMRES coupled method (iter_max=%s)"%max_iter
1461                    p0=self.solve_prec(Bz)
1462               #alfa=(1/self.vol)*util.integrate(util.interpolate(p,escript.Function(self.domain)))
1463                   #p-=alfa
1464                    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)
1465                    #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())
1466                    # solve Au=f-B^*p
1467                    #       A(u-v)=f-B^*p-Av
1468                    #       u=v+(u-v)
1469                p=x[1]
1470            #u=v+self.solve_A(v,p)      
1471            #p=x[1]
1472            #u=x[0]
1473    
1474                  if solver=='PCG':
1475                    #   note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv
1476                    #
1477                    #   with                    Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask)
1478                    #                           A(v-z)= f -Az - B^*p (v-z=0 on fixed_u_mask)
1479                    if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter
1480                    p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p)
1481                u=r[0]  
1482                    print "DIFF=",util.integrate(p)
1483    
1484                  print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)
1485    
1486              return u,p
1487    
1488          def __Msolve(self,r):
1489              return self.solve_prec(r[1])
1490    
1491          def __Msolve2(self,r):
1492              return self.solve_prec(r*1)
1493    
1494          def __Mempty(self,r):
1495              return r
1496    
1497    
1498          def __Aprod(self,p):
1499              # return BA^-1B*p
1500              #solve Av =B^*p as Av =f-Az-B^*(-p)
1501              v=self.solve_A(self.__z,-p)
1502              return ArithmeticTuple(v, self.B(v))
1503    
1504          def __Anext(self,x):
1505              # return next v,p
1506              #solve Av =-B^*p as Av =f-Az-B^*p
1507    
1508          pc=x[1]
1509              v=self.solve_A(self.__z,-pc)
1510          p=self.solve_prec(self.B(v))
1511    
1512              return ArithmeticTuple(v,p)
1513    
1514    
1515          def __Aprod2(self,p):
1516              # return BA^-1B*p
1517              #solve Av =B^*p as Av =f-Az-B^*(-p)
1518          v=self.solve_A(self.__z,-p)
1519              return self.B(v)
1520    
1521          def __Aprod_Newton(self,p):
1522              # return BA^-1B*p - Bz
1523              #solve Av =-B^*p as Av =f-Az-B^*p
1524          v=self.solve_A(self.__z,-p)
1525              return self.B(v-self.__z)
1526    
1527          def __Aprod_Newton2(self,x):
1528              # return BA^-1B*p - Bz
1529              #solve Av =-B^*p as Av =f-Az-B^*p
1530              pc=x[1]
1531          v=self.solve_A(self.__z,-pc)
1532              p=self.solve_prec(self.B(v-self.__z))
1533              return ArithmeticTuple(v,p)
1534    
1535    class SaddlePointProblem(object):
1536       """
1537       This implements a solver for a saddlepoint problem
1538    
1539       M{f(u,p)=0}
1540       M{g(u)=0}
1541    
1542       for u and p. The problem is solved with an inexact Uszawa scheme for p:
1543    
1544       M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})}
1545       M{Q_g (p^{k+1}-p^{k}) =   g(u^{k+1})}
1546    
1547       where Q_f is an approximation of the Jacobiean A_f of f with respect to u  and Q_f is an approximation of
1548       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'
1549       Q_g can be difficult, non-linear conjugate gradient method is applied to solve for p, so Q_g plays
1550       in fact the role of a preconditioner.
1551       """
1552       def __init__(self,verbose=False,*args):
1553           """
1554           initializes the problem
1555    
1556           @param verbose: switches on the printing out some information
1557           @type verbose: C{bool}
1558           @note: this method may be overwritten by a particular saddle point problem
1559           """
1560           if not isinstance(verbose,bool):
1561                raise TypeError("verbose needs to be of type bool.")
1562           self.__verbose=verbose
1563           self.relaxation=1.
1564    
1565       def trace(self,text):
1566           """
1567           prints text if verbose has been set
1568    
1569           @param text: a text message
1570           @type text: C{str}
1571           """
1572           if self.__verbose: print "#s: #s"%(str(self),text)
1573    
1574       def solve_f(self,u,p,tol=1.e-8):
1575           """
1576           solves
1577    
1578           A_f du = f(u,p)
1579    
1580           with tolerance C{tol} and return du. A_f is Jacobiean of f with respect to u.
1581    
1582           @param u: current approximation of u
1583           @type u: L{escript.Data}
1584           @param p: current approximation of p
1585           @type p: L{escript.Data}
1586           @param tol: tolerance expected for du
1587           @type tol: C{float}
1588           @return: increment du
1589           @rtype: L{escript.Data}
1590           @note: this method has to be overwritten by a particular saddle point problem
1591           """
1592           pass
1593    
1594       def solve_g(self,u,tol=1.e-8):
1595           """
1596           solves
1597    
1598           Q_g dp = g(u)
1599    
1600           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.
1601    
1602           @param u: current approximation of u
1603           @type u: L{escript.Data}
1604           @param tol: tolerance expected for dp
1605           @type tol: C{float}
1606           @return: increment dp
1607           @rtype: L{escript.Data}
1608           @note: this method has to be overwritten by a particular saddle point problem
1609           """
1610           pass
1611    
1612       def inner(self,p0,p1):
1613           """
1614           inner product of p0 and p1 approximating p. Typically this returns integrate(p0*p1)
1615           @return: inner product of p0 and p1
1616           @rtype: C{float}
1617           """
1618           pass
1619    
1620       subiter_max=3
1621       def solve(self,u0,p0,tolerance=1.e-6,tolerance_u=None,iter_max=100,accepted_reduction=0.995,relaxation=None):
1622            """
1623            runs the solver
1624    
1625            @param u0: initial guess for C{u}
1626            @type u0: L{esys.escript.Data}
1627            @param p0: initial guess for C{p}
1628            @type p0: L{esys.escript.Data}
1629            @param tolerance: tolerance for relative error in C{u} and C{p}
1630            @type tolerance: positive C{float}
1631            @param tolerance_u: tolerance for relative error in C{u} if different from C{tolerance}
1632            @type tolerance_u: positive C{float}
1633            @param iter_max: maximum number of iteration steps.
1634            @type iter_max: C{int}
1635            @param accepted_reduction: if the norm  g cannot be reduced by C{accepted_reduction} backtracking to adapt the
1636                                       relaxation factor. If C{accepted_reduction=None} no backtracking is used.
1637            @type accepted_reduction: positive C{float} or C{None}
1638            @param relaxation: initial relaxation factor. If C{relaxation==None}, the last relaxation factor is used.
1639            @type relaxation: C{float} or C{None}
1640            """
1641            tol=1.e-2
1642            if tolerance_u==None: tolerance_u=tolerance
1643            if not relaxation==None: self.relaxation=relaxation
1644            if accepted_reduction ==None:
1645                  angle_limit=0.
1646            elif accepted_reduction>=1.:
1647                  angle_limit=0.
1648            else:
1649                  angle_limit=util.sqrt(1-accepted_reduction**2)
1650            self.iter=0
1651            u=u0
1652            p=p0
1653            #
1654            #   initialize things:
1655            #
1656            converged=False
1657            #
1658            #  start loop:
1659            #
1660            #  initial search direction is g
1661            #
1662            while not converged :
1663                if self.iter>iter_max:
1664                    raise ArithmeticError("no convergence after %s steps."%self.iter)
1665                f_new=self.solve_f(u,p,tol)
1666                norm_f_new = util.Lsup(f_new)
1667                u_new=u-f_new
1668                g_new=self.solve_g(u_new,tol)
1669                self.iter+=1
1670                norm_g_new = util.sqrt(self.inner(g_new,g_new))
1671                if norm_f_new==0. and norm_g_new==0.: return u, p
1672                if self.iter>1 and not accepted_reduction==None:
1673                   #
1674                   #   did we manage to reduce the norm of G? I
1675                   #   if not we start a backtracking procedure
1676                   #
1677                   # print "new/old norm = ",norm_g_new, norm_g, norm_g_new/norm_g
1678                   if norm_g_new > accepted_reduction * norm_g:
1679                      sub_iter=0
1680                      s=self.relaxation
1681                      d=g
1682                      g_last=g
1683                      self.trace("    start substepping: f = %s, g = %s, relaxation = %s."%(norm_f_new, norm_g_new, s))
1684                      while sub_iter < self.subiter_max and  norm_g_new > accepted_reduction * norm_g:
1685                         dg= g_new-g_last
1686                         norm_dg=abs(util.sqrt(self.inner(dg,dg))/self.relaxation)
1687                         rad=self.inner(g_new,dg)/self.relaxation
1688                         # print "   ",sub_iter,": rad, norm_dg:",abs(rad), norm_dg*norm_g_new * angle_limit
1689                         # print "   ",sub_iter,": rad, norm_dg:",rad, norm_dg, norm_g_new, norm_g
1690                         if abs(rad) < norm_dg*norm_g_new * angle_limit:
1691                             if sub_iter>0: self.trace("    no further improvements expected from backtracking.")
1692                             break
1693                         r=self.relaxation
1694                         self.relaxation= - rad/norm_dg**2
1695                         s+=self.relaxation
1696                         #####
1697                         # a=g_new+self.relaxation*dg/r
1698                         # print "predicted new norm = ",util.sqrt(self.inner(a,a)),util.sqrt(self.inner(g_new,g_new)), self.relaxation
1699                         #####
1700                         g_last=g_new
1701                         p+=self.relaxation*d
1702                         f_new=self.solve_f(u,p,tol)
1703                         u_new=u-f_new
1704                         g_new=self.solve_g(u_new,tol)
1705                         self.iter+=1
1706                         norm_f_new = util.Lsup(f_new)
1707                         norm_g_new = util.sqrt(self.inner(g_new,g_new))
1708                         # print "   ",sub_iter," new g norm",norm_g_new
1709                         self.trace("    %s th sub-step: f = %s, g = %s, relaxation = %s."%(sub_iter, norm_f_new, norm_g_new, s))
1710                         #
1711                         #   can we expect reduction of g?
1712                         #
1713                         # u_last=u_new
1714                         sub_iter+=1
1715                      self.relaxation=s
1716                #
1717                #  check for convergence:
1718                #
1719                norm_u_new = util.Lsup(u_new)
1720                p_new=p+self.relaxation*g_new
1721                norm_p_new = util.sqrt(self.inner(p_new,p_new))
1722                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))
1723    
1724                if self.iter>1:
1725                   dg2=g_new-g
1726                   df2=f_new-f
1727                   norm_dg2=util.sqrt(self.inner(dg2,dg2))
1728                   norm_df2=util.Lsup(df2)
1729                   # print norm_g_new, norm_g, norm_dg, norm_p, tolerance
1730                   tol_eq_g=tolerance*norm_dg2/(norm_g*abs(self.relaxation))*norm_p_new
1731                   tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new
1732                   if norm_g_new <= tol_eq_g and norm_f_new <= tol_eq_f:
1733                       converged=True
1734                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
1735            self.trace("convergence after %s steps."%self.iter)
1736            return u,p
1737    #   def solve(self,u0,p0,tolerance=1.e-6,iter_max=10,self.relaxation=1.):
1738    #      tol=1.e-2
1739    #      iter=0
1740    #      converged=False
1741    #      u=u0*1.
1742    #      p=p0*1.
1743    #      while not converged and iter<iter_max:
1744    #          du=self.solve_f(u,p,tol)
1745    #          u-=du
1746    #          norm_du=util.Lsup(du)
1747    #          norm_u=util.Lsup(u)
1748    #        
1749    #          dp=self.relaxation*self.solve_g(u,tol)
1750    #          p+=dp
1751    #          norm_dp=util.sqrt(self.inner(dp,dp))
1752    #          norm_p=util.sqrt(self.inner(p,p))
1753    #          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)
1754    #          iter+=1
1755    #
1756    #          converged = (norm_du <= tolerance*norm_u) and  (norm_dp <= tolerance*norm_p)
1757    #      if converged:
1758    #          print "convergence after %s steps."%iter
1759    #      else:
1760    #          raise ArithmeticError("no convergence after %s steps."%iter)
1761    #
1762    #      return u,p
1763              
1764    def MaskFromBoundaryTag(function_space,*tags):
1765       """
1766       create a mask on the given function space which one for samples
1767       that touch the boundary tagged by tags.
1768    
1769       usage: m=MaskFromBoundaryTag(Solution(domain),"left", "right")
1770    
1771       @param function_space: a given function space
1772       @type function_space: L{escript.FunctionSpace}
1773       @param tags: boundray tags
1774       @type tags: C{str}
1775       @return: a mask which marks samples used by C{function_space} that are touching the
1776                boundary tagged by any of the given tags.
1777       @rtype: L{escript.Data} of rank 0
1778       """
1779       pde=linearPDEs.LinearPDE(function_space.getDomain(),numEquations=1, numSolutions=1)
1780       d=escript.Scalar(0.,escript.FunctionOnBoundary(function_space.getDomain()))
1781       for t in tags: d.setTaggedValue(t,1.)
1782       pde.setValue(y=d)
1783       out=util.whereNonZero(pde.getRightHandSide())
1784       if out.getFunctionSpace() == function_space:
1785          return out
1786       else:
1787          return util.whereNonZero(util.interpolate(out,function_space))
1788    
1789    
1790    

Legend:
Removed from v.396  
changed lines
  Added in v.1552

  ViewVC Help
Powered by ViewVC 1.1.26