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

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

  ViewVC Help
Powered by ViewVC 1.1.26