/[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

trunk/escript/py_src/pdetools.py revision 1122 by gross, Tue May 1 03:21:04 2007 UTC temp/escript/py_src/pdetools.py revision 1387 by trankine, Fri Jan 11 07:45:26 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 32  import escript Line 46  import escript
46  import linearPDEs  import linearPDEs
47  import numarray  import numarray
48  import util  import util
49    import math
50    
51  class TimeIntegrationManager:  class TimeIntegrationManager:
52    """    """
# Line 393  class Locator: Line 408  class Locator:
408          else:          else:
409             return data             return data
410    
411    class SolverSchemeException(Exception):
412       """
413       exceptions thrown by solvers
414       """
415       pass
416    
417    class IndefinitePreconditioner(SolverSchemeException):
418       """
419       the preconditioner is not positive definite.
420       """
421       pass
422    class MaxIterReached(SolverSchemeException):
423       """
424       maxium number of iteration steps is reached.
425       """
426       pass
427    class IterationBreakDown(SolverSchemeException):
428       """
429       iteration scheme econouters an incurable breakdown.
430       """
431       pass
432    class NegativeNorm(SolverSchemeException):
433       """
434       a norm calculation returns a negative norm.
435       """
436       pass
437    
438    class IterationHistory(object):
439       """
440       The IterationHistory class is used to define a stopping criterium. It keeps track of the
441       residual norms. The stoppingcriterium indicates termination if the residual norm has been reduced by
442       a given tolerance.
443       """
444       def __init__(self,tolerance=math.sqrt(util.EPSILON),verbose=False):
445          """
446          Initialization
447    
448          @param tolerance: tolerance
449          @type tolerance: positive C{float}
450          @param verbose: switches on the printing out some information
451          @type verbose: C{bool}
452          """
453          if not tolerance>0.:
454              raise ValueError,"tolerance needs to be positive."
455          self.tolerance=tolerance
456          self.verbose=verbose
457          self.history=[]
458       def stoppingcriterium(self,norm_r,r,x):
459           """
460           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.
461    
462          
463           @param norm_r: current residual norm
464           @type norm_r: non-negative C{float}
465           @param r: current residual (not used)
466           @param x: current solution approximation (not used)
467           @return: C{True} is the stopping criterium is fullfilled. Otherwise C{False} is returned.
468           @rtype: C{bool}
469    
470           """
471           self.history.append(norm_r)
472           if self.verbose: print "iter: %s:  inner(rhat,r) = %e"%(len(self.history)-1, self.history[-1])
473           return self.history[-1]<=self.tolerance * self.history[0]
474    
475    def PCG(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
476       """
477       Solver for
478    
479       M{Ax=b}
480    
481       with a symmetric and positive definite operator A (more details required!).
482       It uses the conjugate gradient method with preconditioner M providing an approximation of A.
483    
484       The iteration is terminated if the C{stoppingcriterium} function return C{True}.
485    
486       For details on the preconditioned conjugate gradient method see the book:
487    
488       Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
489       T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
490       C. Romine, and H. van der Vorst.
491    
492       @param b: the right hand side of the liner system. C{b} is altered.
493       @type b: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
494       @param Aprod: returns the value Ax
495       @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}.
496       @param Msolve: solves Mx=r
497       @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
498    type like argument C{x}.
499       @param bilinearform: inner product C{<x,r>}
500       @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}.
501       @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.
502       @type stoppingcriterium: function that returns C{True} or C{False}
503       @param x: an initial guess for the solution. If no C{x} is given 0*b is used.
504       @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
505       @param iter_max: maximum number of iteration steps.
506       @type iter_max: C{int}
507       @return: the solution approximation and the corresponding residual
508       @rtype: C{tuple}
509       @warning: C{b} and C{x} are altered.
510       """
511       iter=0
512       if x==None:
513          x=0*b
514       else:
515          b += (-1)*Aprod(x)
516       r=b
517       rhat=Msolve(r)
518       d = rhat
519       rhat_dot_r = bilinearform(rhat, r)
520       if rhat_dot_r<0: raise NegativeNorm,"negative norm."
521    
522       while not stoppingcriterium(math.sqrt(rhat_dot_r),r,x):
523           iter+=1
524           if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
525    
526           q=Aprod(d)
527           alpha = rhat_dot_r / bilinearform(d, q)
528           x += alpha * d
529           r += (-alpha) * q
530    
531           rhat=Msolve(r)
532           rhat_dot_r_new = bilinearform(rhat, r)
533           beta = rhat_dot_r_new / rhat_dot_r
534           rhat+=beta * d
535           d=rhat
536    
537           rhat_dot_r = rhat_dot_r_new
538           if rhat_dot_r<0: raise NegativeNorm,"negative norm."
539    
540       return x,r
541    
542    class ArithmeticTuple(object):
543       """
544       tuple supporting inplace update x+=y and scaling x=a*y where x,y is an ArithmeticTuple and a is a float.
545    
546       example of usage:
547    
548       from esys.escript import Data
549       from numarray import array
550       a=Data(...)
551       b=array([1.,4.])
552       x=ArithmeticTuple(a,b)
553       y=5.*x
554    
555       """
556       def __init__(self,*args):
557           """
558           initialize object with elements args.
559    
560           @param args: tuple of object that support implace add (x+=y) and scaling (x=a*y)
561           """
562           self.__items=list(args)
563    
564       def __len__(self):
565           """
566           number of items
567    
568           @return: number of items
569           @rtype: C{int}
570           """
571           return len(self.__items)
572    
573       def __getitem__(self,index):
574           """
575           get an item
576    
577           @param index: item to be returned
578           @type index: C{int}
579           @return: item with index C{index}
580           """
581           return self.__items.__getitem__(index)
582    
583       def __mul__(self,other):
584           """
585           scaling from the right
586    
587           @param other: scaling factor
588           @type other: C{float}
589           @return: itemwise self*other
590           @rtype: L{ArithmeticTuple}
591           """
592           out=[]
593           for i in range(len(self)):
594               out.append(self[i]*other)
595           return ArithmeticTuple(*tuple(out))
596    
597       def __rmul__(self,other):
598           """
599           scaling from the left
600    
601           @param other: scaling factor
602           @type other: C{float}
603           @return: itemwise other*self
604           @rtype: L{ArithmeticTuple}
605           """
606           out=[]
607           for i in range(len(self)):
608               out.append(other*self[i])
609           return ArithmeticTuple(*tuple(out))
610    
611       def __iadd__(self,other):
612           """
613           in-place add of other to self
614    
615           @param other: increment
616           @type other: C{ArithmeticTuple}
617           """
618           if len(self) != len(other):
619               raise ValueError,"tuple length must match."
620           for i in range(len(self)):
621               self.__items[i]+=other[i]
622           return self
623    
624  class SaddlePointProblem(object):  class SaddlePointProblem(object):
625     """     """
626     This implements a solver for a saddlepoint problem     This implements a solver for a saddlepoint problem
# Line 580  class SaddlePointProblem(object): Line 808  class SaddlePointProblem(object):
808              norm_u_new = util.Lsup(u_new)              norm_u_new = util.Lsup(u_new)
809              p_new=p+self.relaxation*g_new              p_new=p+self.relaxation*g_new
810              norm_p_new = util.sqrt(self.inner(p_new,p_new))              norm_p_new = util.sqrt(self.inner(p_new,p_new))
811              self.trace("%s th step: f/u = %s, g/p = %s, relaxation = %s."%(self.iter,norm_f_new/norm_u_new, norm_g_new/norm_p_new, self.relaxation))              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))
812    
813              if self.iter>1:              if self.iter>1:
814                 dg2=g_new-g                 dg2=g_new-g
# Line 622  class SaddlePointProblem(object): Line 850  class SaddlePointProblem(object):
850  #  #
851  #      return u,p  #      return u,p
852                        
853  # vim: expandtab shiftwidth=4:  def MaskFromBoundaryTag(function_space,*tags):
854       """
855       create a mask on the given function space which one for samples
856       that touch the boundary tagged by tags.
857    
858       usage: m=MaskFromBoundaryTag(Solution(domain),"left", "right")
859    
860       @param function_space: a given function space
861       @type function_space: L{escript.FunctionSpace}
862       @param tags: boundray tags
863       @type tags: C{str}
864       @return: a mask which marks samples used by C{function_space} that are touching the
865                boundary tagged by any of the given tags.
866       @rtype: L{escript.Data} of rank 0
867       """
868       pde=linearPDEs.LinearPDE(function_space.getDomain(),numEquations=1, numSolutions=1)
869       d=escript.Scalar(0.,escript.FunctionOnBoundary(function_space.getDomain()))
870       for t in tags: d.setTaggedValue(t,1.)
871       pde.setValue(y=d)
872       out=util.whereNonZero(pde.getRightHandSide())
873       if out.getFunctionSpace() == function_space:
874          return out
875       else:
876          return util.whereNonZero(util.interpolate(out,function_space))
877    

Legend:
Removed from v.1122  
changed lines
  Added in v.1387

  ViewVC Help
Powered by ViewVC 1.1.26