/[escript]/trunk/escript/py_src/pdetools.py
ViewVC logotype

Diff of /trunk/escript/py_src/pdetools.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 351 by gross, Tue Dec 13 09:12:15 2005 UTC revision 1465 by artak, Wed Apr 2 03:28:25 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,Float32,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 49  class TimeIntegrationManager: Line 88  class TimeIntegrationManager:
88    
89    def getTime(self):    def getTime(self):
90        return self.__t        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):    def checkin(self,dt,*values):
99        """        """
# Line 84  class TimeIntegrationManager: Line 129  class TimeIntegrationManager:
129           return out[0]           return out[0]
130        else:        else:
131           return out           return out
132    
133    
134  class Projector:  class Projector:
135    """    """
# Line 106  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 143  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 159  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 178  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 204  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 PCG(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
481       """
482       Solver for
483    
484       M{Ax=b}
485    
486       with a symmetric and positive definite operator A (more details required!).
487       It uses the conjugate gradient method with preconditioner M providing an approximation of A.
488    
489       The iteration is terminated if the C{stoppingcriterium} function return C{True}.
490    
491       For details on the preconditioned conjugate gradient method see the book:
492    
493       Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
494       T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
495       C. Romine, and H. van der Vorst.
496    
497       @param b: the right hand side of the liner system. C{b} is altered.
498       @type b: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
499       @param Aprod: returns the value Ax
500       @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}.
501       @param Msolve: solves Mx=r
502       @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
503    type like argument C{x}.
504       @param bilinearform: inner product C{<x,r>}
505       @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}.
506       @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.
507       @type stoppingcriterium: function that returns C{True} or C{False}
508       @param x: an initial guess for the solution. If no C{x} is given 0*b is used.
509       @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
510       @param iter_max: maximum number of iteration steps.
511       @type iter_max: C{int}
512       @return: the solution approximation and the corresponding residual
513       @rtype: C{tuple}
514       @warning: C{b} and C{x} are altered.
515       """
516       iter=0
517       if x==None:
518          x=0*b
519       else:
520          b += (-1)*Aprod(x)
521       r=b
522       rhat=Msolve(r)
523       d = rhat
524       rhat_dot_r = bilinearform(rhat, r)
525       if rhat_dot_r<0: raise NegativeNorm,"negative norm."
526    
527       while not stoppingcriterium(math.sqrt(rhat_dot_r),r,x):
528           iter+=1
529           if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
530    
531           q=Aprod(d)
532           alpha = rhat_dot_r / bilinearform(d, q)
533           x += alpha * d
534           r += (-alpha) * q
535    
536           rhat=Msolve(r)
537           rhat_dot_r_new = bilinearform(rhat, r)
538           beta = rhat_dot_r_new / rhat_dot_r
539           rhat+=beta * d
540           d=rhat
541    
542           rhat_dot_r = rhat_dot_r_new
543           if rhat_dot_r<0: raise NegativeNorm,"negative norm."
544    
545       return x,r
546    
547    
548    ############################
549    # Added by Artak
550    #################################3
551    
552    #Apply a sequence of k Givens rotations, used within gmres codes
553    # vrot=givapp(c, s, vin, k)
554    def givapp(c,s,vin,k):
555        vrot=vin
556        for i in range(k+1):
557            w1=c[i]*vrot[i]-s[i]*vrot[i+1]
558        w2=s[i]*vrot[i]+c[i]*vrot[i+1]
559            vrot[i:i+2]=w1,w2
560        return vrot
561    
562    def GMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
563    
564       from numarray import dot
565      
566       v0=b
567       iter=0
568       if x==None:
569          x=0*b
570       else:
571          b += (-1.)*Aprod(x)
572       r=b
573    
574       rhat=Msolve(r)
575      
576       rhat_dot_r = bilinearform(rhat, rhat)
577       norm_r=math.sqrt(rhat_dot_r)
578    
579       if rhat_dot_r<0: raise NegativeNorm,"negative norm."
580      
581       h=zeros((iter_max,iter_max),Float32)
582       c=zeros(iter_max,Float32)
583       s=zeros(iter_max,Float32)
584       g=zeros(iter_max,Float32)
585       v=[]
586    
587       v.append(rhat/norm_r)
588       rho=norm_r
589       g[0]=rho
590    
591       while not stoppingcriterium(rho,rho,norm_r):
592    
593        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
594    
595        
596        vhat=Aprod(v[iter])
597        p=Msolve(vhat)  
598    
599        v.append(p)
600    
601        v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
602    
603    # Modified Gram-Schmidt
604        for j in range(iter+1):
605          h[j][iter]=bilinearform(v[j],v[iter+1])  
606          v[iter+1]+=(-1.)*h[j][iter]*v[j]
607          
608        h[iter+1][iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
609        v_norm2=h[iter+1][iter]
610    
611    
612    # Reorthogonalize if needed
613        if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
614         for j in range(iter+1):
615            hr=bilinearform(v[j],v[iter+1])
616                h[j][iter]=h[j][iter]+hr #vhat
617                v[iter+1] +=(-1.)*hr*v[j]
618    
619         v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
620         h[iter+1][iter]=v_norm2
621    
622    #   watch out for happy breakdown
623            if v_norm2 != 0:
624             v[iter+1]=v[iter+1]/h[iter+1][iter]
625    
626    #   Form and store the information for the new Givens rotation
627        if iter > 0 :
628            hhat=[]
629            for i in range(iter+1) : hhat.append(h[i][iter])
630            hhat=givapp(c[0:iter],s[0:iter],hhat,iter-1);
631                for i in range(iter+1) : h[i][iter]=hhat[i]
632    
633        mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])
634        if mu!=0 :
635            c[iter]=h[iter][iter]/mu
636            s[iter]=-h[iter+1][iter]/mu
637            h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter]
638            h[iter+1][iter]=0.0
639            g[iter:iter+2]=givapp(c[iter],s[iter],g[iter:iter+2],0)
640    
641    # Update the residual norm
642            rho=abs(g[iter+1])
643        iter+=1
644    
645    # At this point either iter > iter_max or rho < tol.
646    # It's time to compute x and leave.        
647    
648       if iter > 0 :
649         y=zeros(iter,Float32)  
650         y[iter-1] = g[iter-1] / h[iter-1][iter-1]
651         if iter > 1 :  
652            i=iter-2  
653            while i>=0 :
654              y[i] = ( g[i] - dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]
655              i=i-1
656         xhat=v[iter-1]*y[iter-1]
657         for i in range(iter-1):
658        xhat += v[i]*y[i]
659       else : xhat=v[0]
660        
661       x += xhat
662    
663       return x
664        
665    #############################################
666    
667    class ArithmeticTuple(object):
668       """
669       tuple supporting inplace update x+=y and scaling x=a*y where x,y is an ArithmeticTuple and a is a float.
670    
671       example of usage:
672    
673       from esys.escript import Data
674       from numarray import array
675       a=Data(...)
676       b=array([1.,4.])
677       x=ArithmeticTuple(a,b)
678       y=5.*x
679    
680       """
681       def __init__(self,*args):
682           """
683           initialize object with elements args.
684    
685           @param args: tuple of object that support implace add (x+=y) and scaling (x=a*y)
686           """
687           self.__items=list(args)
688    
689       def __len__(self):
690           """
691           number of items
692    
693           @return: number of items
694           @rtype: C{int}
695           """
696           return len(self.__items)
697    
698       def __getitem__(self,index):
699           """
700           get an item
701    
702           @param index: item to be returned
703           @type index: C{int}
704           @return: item with index C{index}
705           """
706           return self.__items.__getitem__(index)
707    
708       def __mul__(self,other):
709           """
710           scaling from the right
711    
712           @param other: scaling factor
713           @type other: C{float}
714           @return: itemwise self*other
715           @rtype: L{ArithmeticTuple}
716           """
717           out=[]
718           for i in range(len(self)):
719               out.append(self[i]*other)
720           return ArithmeticTuple(*tuple(out))
721    
722       def __rmul__(self,other):
723           """
724           scaling from the left
725    
726           @param other: scaling factor
727           @type other: C{float}
728           @return: itemwise other*self
729           @rtype: L{ArithmeticTuple}
730           """
731           out=[]
732           for i in range(len(self)):
733               out.append(other*self[i])
734           return ArithmeticTuple(*tuple(out))
735    
736    #########################
737    # Added by Artak
738    #########################
739       def __div__(self,other):
740           """
741           dividing from the right
742    
743           @param other: scaling factor
744           @type other: C{float}
745           @return: itemwise self/other
746           @rtype: L{ArithmeticTuple}
747           """
748           out=[]
749           for i in range(len(self)):
750               out.append(self[i]/other)
751           return ArithmeticTuple(*tuple(out))
752    
753       def __rdiv__(self,other):
754           """
755           dividing from the left
756    
757           @param other: scaling factor
758           @type other: C{float}
759           @return: itemwise other/self
760           @rtype: L{ArithmeticTuple}
761           """
762           out=[]
763           for i in range(len(self)):
764               out.append(other/self[i])
765           return ArithmeticTuple(*tuple(out))
766      
767    ##########################################33
768    
769       def __iadd__(self,other):
770           """
771           in-place add of other to self
772    
773           @param other: increment
774           @type other: C{ArithmeticTuple}
775           """
776           if len(self) != len(other):
777               raise ValueError,"tuple length must match."
778           for i in range(len(self)):
779               self.__items[i]+=other[i]
780           return self
781    
782    class HomogeneousSaddlePointProblem(object):
783          """
784          This provides a framwork for solving homogeneous saddle point problem of the form
785    
786                 Av+B^*p=f
787                 Bv    =0
788    
789          for the unknowns v and p and given operators A and B and given right hand side f.
790          B^* is the adjoint operator of B is the given inner product.
791    
792          """
793          def __init__(self,**kwargs):
794            self.setTolerance()
795            self.setToleranceReductionFactor()
796    
797          def initialize(self):
798            """
799            initialize the problem (overwrite)
800            """
801            pass
802          def B(self,v):
803             """
804             returns Bv (overwrite)
805             @rtype: equal to the type of p
806    
807             @note: boundary conditions on p should be zero!
808             """
809             pass
810    
811          def inner(self,p0,p1):
812             """
813             returns inner product of two element p0 and p1  (overwrite)
814            
815             @type p0: equal to the type of p
816             @type p1: equal to the type of p
817             @rtype: C{float}
818    
819             @rtype: equal to the type of p
820             """
821             pass
822    
823          def solve_A(self,u,p):
824             """
825             solves Av=f-Au-B^*p with accuracy self.getReducedTolerance() (overwrite)
826    
827             @rtype: equal to the type of v
828             @note: boundary conditions on v should be zero!
829             """
830             pass
831    
832          def solve_prec(self,p):
833             """
834             provides a preconditioner for BA^{-1}B^* with accuracy self.getReducedTolerance() (overwrite)
835    
836             @rtype: equal to the type of p
837             """
838             pass
839    
840          def stoppingcriterium(self,Bv,v,p):
841             """
842             returns a True if iteration is terminated. (overwrite)
843    
844             @rtype: C{bool}
845             """
846             pass
847                
848          def __inner(self,p,r):
849             return self.inner(p,r[1])
850    
851          def __inner_p(self,p1,p2):
852             return self.inner(p1,p2)
853    
854          def __stoppingcriterium(self,norm_r,r,p):
855              return self.stoppingcriterium(r[1],r[0],p)
856    
857          def __stoppingcriterium_GMRES(self,norm_r,rho,r):
858              return self.stoppingcriterium_GMRES(rho,r)
859    
860          def setTolerance(self,tolerance=1.e-8):
861                  self.__tol=tolerance
862          def getTolerance(self):
863                  return self.__tol
864          def setToleranceReductionFactor(self,reduction=0.01):
865                  self.__reduction=reduction
866          def getSubProblemTolerance(self):
867                  return self.__reduction*self.getTolerance()
868    
869          def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='GMRES'):
870                  """
871                  solves the saddle point problem using initial guesses v and p.
872    
873                  @param max_iter: maximum number of iteration steps.
874                  """
875                  self.verbose=verbose
876                  self.show_details=show_details and self.verbose
877    
878              # Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask)      
879    
880              
881              self.__z=v+self.solve_A(v,p*0)
882    
883                  Bz=self.B(self.__z)
884                  #
885              #   solve BA^-1B^*p = Bz
886                  #
887                  #   note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv
888                  #
889                  #   with                    Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask)
890                  #                           A(v-z)=Az-B^*p-Az = f -Az - B^*p (v-z=0 on fixed_u_mask)
891                  #
892                  self.iter=0
893              if solver=='GMRES':      
894                    if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter
895                    p=GMRES(Bz,self.__Aprod_GMRES,self.__Msolve_GMRES,self.__inner_p,self.__stoppingcriterium_GMRES,iter_max=max_iter, x=p*1)
896            u=v+self.solve_A(v,p)
897        
898                  else:
899                    if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter
900                    p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p*1)
901                u=r[0]  
902                    u=v+self.solve_A(v,p)     # Lutz to check !!!!!
903    
904              return u,p
905    
906          def __Msolve(self,r):
907              return self.solve_prec(r[1])
908    
909          def __Msolve_GMRES(self,r):
910              return self.solve_prec(r)
911    
912    
913          def __Aprod(self,p):
914              # return BA^-1B*p
915              #solve Av =-B^*p as Av =f-Az-B^*p
916              v=self.solve_A(self.__z,p)
917              return ArithmeticTuple(v, self.B(v))
918    
919          def __Aprod_GMRES(self,p):
920              # return BA^-1B*p
921              #solve Av =-B^*p as Av =f-Az-B^*p
922          v=self.solve_A(self.__z,p)
923              return self.B(v)
924    
925    class SaddlePointProblem(object):
926       """
927       This implements a solver for a saddlepoint problem
928    
929       M{f(u,p)=0}
930       M{g(u)=0}
931    
932       for u and p. The problem is solved with an inexact Uszawa scheme for p:
933    
934       M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})}
935       M{Q_g (p^{k+1}-p^{k}) =   g(u^{k+1})}
936    
937       where Q_f is an approximation of the Jacobiean A_f of f with respect to u  and Q_f is an approximation of
938       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'
939       Q_g can be difficult, non-linear conjugate gradient method is applied to solve for p, so Q_g plays
940       in fact the role of a preconditioner.
941       """
942       def __init__(self,verbose=False,*args):
943           """
944           initializes the problem
945    
946           @param verbose: switches on the printing out some information
947           @type verbose: C{bool}
948           @note: this method may be overwritten by a particular saddle point problem
949           """
950           if not isinstance(verbose,bool):
951                raise TypeError("verbose needs to be of type bool.")
952           self.__verbose=verbose
953           self.relaxation=1.
954    
955       def trace(self,text):
956           """
957           prints text if verbose has been set
958    
959           @param text: a text message
960           @type text: C{str}
961           """
962           if self.__verbose: print "%s: %s"%(str(self),text)
963    
964       def solve_f(self,u,p,tol=1.e-8):
965           """
966           solves
967    
968           A_f du = f(u,p)
969    
970           with tolerance C{tol} and return du. A_f is Jacobiean of f with respect to u.
971    
972           @param u: current approximation of u
973           @type u: L{escript.Data}
974           @param p: current approximation of p
975           @type p: L{escript.Data}
976           @param tol: tolerance expected for du
977           @type tol: C{float}
978           @return: increment du
979           @rtype: L{escript.Data}
980           @note: this method has to be overwritten by a particular saddle point problem
981           """
982           pass
983    
984       def solve_g(self,u,tol=1.e-8):
985           """
986           solves
987    
988           Q_g dp = g(u)
989    
990           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.
991    
992           @param u: current approximation of u
993           @type u: L{escript.Data}
994           @param tol: tolerance expected for dp
995           @type tol: C{float}
996           @return: increment dp
997           @rtype: L{escript.Data}
998           @note: this method has to be overwritten by a particular saddle point problem
999           """
1000           pass
1001    
1002       def inner(self,p0,p1):
1003           """
1004           inner product of p0 and p1 approximating p. Typically this returns integrate(p0*p1)
1005           @return: inner product of p0 and p1
1006           @rtype: C{float}
1007           """
1008           pass
1009    
1010       subiter_max=3
1011       def solve(self,u0,p0,tolerance=1.e-6,tolerance_u=None,iter_max=100,accepted_reduction=0.995,relaxation=None):
1012            """
1013            runs the solver
1014    
1015            @param u0: initial guess for C{u}
1016            @type u0: L{esys.escript.Data}
1017            @param p0: initial guess for C{p}
1018            @type p0: L{esys.escript.Data}
1019            @param tolerance: tolerance for relative error in C{u} and C{p}
1020            @type tolerance: positive C{float}
1021            @param tolerance_u: tolerance for relative error in C{u} if different from C{tolerance}
1022            @type tolerance_u: positive C{float}
1023            @param iter_max: maximum number of iteration steps.
1024            @type iter_max: C{int}
1025            @param accepted_reduction: if the norm  g cannot be reduced by C{accepted_reduction} backtracking to adapt the
1026                                       relaxation factor. If C{accepted_reduction=None} no backtracking is used.
1027            @type accepted_reduction: positive C{float} or C{None}
1028            @param relaxation: initial relaxation factor. If C{relaxation==None}, the last relaxation factor is used.
1029            @type relaxation: C{float} or C{None}
1030            """
1031            tol=1.e-2
1032            if tolerance_u==None: tolerance_u=tolerance
1033            if not relaxation==None: self.relaxation=relaxation
1034            if accepted_reduction ==None:
1035                  angle_limit=0.
1036            elif accepted_reduction>=1.:
1037                  angle_limit=0.
1038            else:
1039                  angle_limit=util.sqrt(1-accepted_reduction**2)
1040            self.iter=0
1041            u=u0
1042            p=p0
1043            #
1044            #   initialize things:
1045            #
1046            converged=False
1047            #
1048            #  start loop:
1049            #
1050            #  initial search direction is g
1051            #
1052            while not converged :
1053                if self.iter>iter_max:
1054                    raise ArithmeticError("no convergence after %s steps."%self.iter)
1055                f_new=self.solve_f(u,p,tol)
1056                norm_f_new = util.Lsup(f_new)
1057                u_new=u-f_new
1058                g_new=self.solve_g(u_new,tol)
1059                self.iter+=1
1060                norm_g_new = util.sqrt(self.inner(g_new,g_new))
1061                if norm_f_new==0. and norm_g_new==0.: return u, p
1062                if self.iter>1 and not accepted_reduction==None:
1063                   #
1064                   #   did we manage to reduce the norm of G? I
1065                   #   if not we start a backtracking procedure
1066                   #
1067                   # print "new/old norm = ",norm_g_new, norm_g, norm_g_new/norm_g
1068                   if norm_g_new > accepted_reduction * norm_g:
1069                      sub_iter=0
1070                      s=self.relaxation
1071                      d=g
1072                      g_last=g
1073                      self.trace("    start substepping: f = %s, g = %s, relaxation = %s."%(norm_f_new, norm_g_new, s))
1074                      while sub_iter < self.subiter_max and  norm_g_new > accepted_reduction * norm_g:
1075                         dg= g_new-g_last
1076                         norm_dg=abs(util.sqrt(self.inner(dg,dg))/self.relaxation)
1077                         rad=self.inner(g_new,dg)/self.relaxation
1078                         # print "   ",sub_iter,": rad, norm_dg:",abs(rad), norm_dg*norm_g_new * angle_limit
1079                         # print "   ",sub_iter,": rad, norm_dg:",rad, norm_dg, norm_g_new, norm_g
1080                         if abs(rad) < norm_dg*norm_g_new * angle_limit:
1081                             if sub_iter>0: self.trace("    no further improvements expected from backtracking.")
1082                             break
1083                         r=self.relaxation
1084                         self.relaxation= - rad/norm_dg**2
1085                         s+=self.relaxation
1086                         #####
1087                         # a=g_new+self.relaxation*dg/r
1088                         # print "predicted new norm = ",util.sqrt(self.inner(a,a)),util.sqrt(self.inner(g_new,g_new)), self.relaxation
1089                         #####
1090                         g_last=g_new
1091                         p+=self.relaxation*d
1092                         f_new=self.solve_f(u,p,tol)
1093                         u_new=u-f_new
1094                         g_new=self.solve_g(u_new,tol)
1095                         self.iter+=1
1096                         norm_f_new = util.Lsup(f_new)
1097                         norm_g_new = util.sqrt(self.inner(g_new,g_new))
1098                         # print "   ",sub_iter," new g norm",norm_g_new
1099                         self.trace("    %s th sub-step: f = %s, g = %s, relaxation = %s."%(sub_iter, norm_f_new, norm_g_new, s))
1100                         #
1101                         #   can we expect reduction of g?
1102                         #
1103                         # u_last=u_new
1104                         sub_iter+=1
1105                      self.relaxation=s
1106                #
1107                #  check for convergence:
1108                #
1109                norm_u_new = util.Lsup(u_new)
1110                p_new=p+self.relaxation*g_new
1111                norm_p_new = util.sqrt(self.inner(p_new,p_new))
1112                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))
1113    
1114                if self.iter>1:
1115                   dg2=g_new-g
1116                   df2=f_new-f
1117                   norm_dg2=util.sqrt(self.inner(dg2,dg2))
1118                   norm_df2=util.Lsup(df2)
1119                   # print norm_g_new, norm_g, norm_dg, norm_p, tolerance
1120                   tol_eq_g=tolerance*norm_dg2/(norm_g*abs(self.relaxation))*norm_p_new
1121                   tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new
1122                   if norm_g_new <= tol_eq_g and norm_f_new <= tol_eq_f:
1123                       converged=True
1124                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
1125            self.trace("convergence after %s steps."%self.iter)
1126            return u,p
1127    #   def solve(self,u0,p0,tolerance=1.e-6,iter_max=10,self.relaxation=1.):
1128    #      tol=1.e-2
1129    #      iter=0
1130    #      converged=False
1131    #      u=u0*1.
1132    #      p=p0*1.
1133    #      while not converged and iter<iter_max:
1134    #          du=self.solve_f(u,p,tol)
1135    #          u-=du
1136    #          norm_du=util.Lsup(du)
1137    #          norm_u=util.Lsup(u)
1138    #        
1139    #          dp=self.relaxation*self.solve_g(u,tol)
1140    #          p+=dp
1141    #          norm_dp=util.sqrt(self.inner(dp,dp))
1142    #          norm_p=util.sqrt(self.inner(p,p))
1143    #          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)
1144    #          iter+=1
1145    #
1146    #          converged = (norm_du <= tolerance*norm_u) and  (norm_dp <= tolerance*norm_p)
1147    #      if converged:
1148    #          print "convergence after %s steps."%iter
1149    #      else:
1150    #          raise ArithmeticError("no convergence after %s steps."%iter)
1151    #
1152    #      return u,p
1153              
1154    def MaskFromBoundaryTag(function_space,*tags):
1155       """
1156       create a mask on the given function space which one for samples
1157       that touch the boundary tagged by tags.
1158    
1159       usage: m=MaskFromBoundaryTag(Solution(domain),"left", "right")
1160    
1161       @param function_space: a given function space
1162       @type function_space: L{escript.FunctionSpace}
1163       @param tags: boundray tags
1164       @type tags: C{str}
1165       @return: a mask which marks samples used by C{function_space} that are touching the
1166                boundary tagged by any of the given tags.
1167       @rtype: L{escript.Data} of rank 0
1168       """
1169       pde=linearPDEs.LinearPDE(function_space.getDomain(),numEquations=1, numSolutions=1)
1170       d=escript.Scalar(0.,escript.FunctionOnBoundary(function_space.getDomain()))
1171       for t in tags: d.setTaggedValue(t,1.)
1172       pde.setValue(y=d)
1173       out=util.whereNonZero(pde.getRightHandSide())
1174       if out.getFunctionSpace() == function_space:
1175          return out
1176       else:
1177          return util.whereNonZero(util.interpolate(out,function_space))
1178    
1179    
1180    

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

  ViewVC Help
Powered by ViewVC 1.1.26