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

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

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

revision 1137 by gross, Thu May 10 08:11:31 2007 UTC revision 1312 by ksteube, Mon Sep 24 06:18:44 2007 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 130  class Projector: Line 145  class Projector:
145      self.__pde.setSymmetryOn()      self.__pde.setSymmetryOn()
146      self.__pde.setReducedOrderTo(reduce)      self.__pde.setReducedOrderTo(reduce)
147      self.__pde.setValue(D = 1.)      self.__pde.setValue(D = 1.)
148        return
149    
150    def __call__(self, input_data):    def __call__(self, input_data):
151      """      """
# Line 392  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    def PCG(b,x,Aprod,Msolve,bilinearform, norm, verbose=True, iter_max=100, tolerance=math.sqrt(util.EPSILON)):
439       """
440       Solver for
441    
442       M{Ax=b}
443    
444       with a symmetric and positive definite operator A (more details required!).
445       It uses the conjugate gradient method with preconditioner M providing an approximation of A.
446    
447       The iteration is terminated if
448    
449       M{norm(r) <= tolerance * norm(b)}
450    
451       where C{norm()} defines a norm and
452    
453       M{r = b-Ax}
454    
455       the residual.
456    
457       For details on the preconditioned conjugate gradient method see the book:
458    
459       Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
460       T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
461       C. Romine, and H. van der Vorst.
462    
463       @param b: the right hand side of the liner system. C{b} is altered.
464       @type b: any object type R supporting inplace add (x+=y) and scaling (x=scalar*y)
465       @param x: an initial guess for the solution
466       @type x: any object type S supporting inplace add (x+=y), scaling (x=scalar*y)
467       @param Aprod: returns the value Ax
468       @type Aprod: function C{Aprod(x)} where C{x} is of object type S. The returned object needs to be of type R.
469       @param Msolve: solves Mx=r
470       @type Msolve: function C{Msolve(r)} where C{r} is of object type R. The returned object needs to be of tupe S.
471       @param bilinearform: inner product C{<x,r>}
472       @type bilinearform: function C{bilinearform(x,r)} where C{x} is of object type S and C{r} is of object type R. The returned value is a C{float}.
473       @param norm: norm calculation for the residual C{r=b-Ax}.
474       @type norm: function C{norm(r)} where C{r} is of object type R. The returned value is a C{float}.
475       @param verbose: switches on the printing out some information
476       @type verbose: C{bool}
477       @param iter_max: maximum number of iteration steps.
478       @type iter_max: C{int}
479       @param tolerance: tolerance
480       @type tolerance: positive C{float}
481       @return: the solution apprximation and the corresponding residual
482       @rtype: C{tuple} of an S type and and an R type object.A
483       @warning: C{b} ans C{x} are altered.
484       """
485       if verbose:
486            print "Enter PCG for solving Ax=b\n\t iter_max =%s\t tolerance   =%e"%(iter_max,tolerance)
487       iter=0
488       normb = norm(b)
489       if normb<0: raise NegativeNorm
490    
491       b += (-1)*Aprod(x)
492       r=b
493       rhat=Msolve(r)
494       d = rhat;
495       rhat_dot_r = bilinearform(rhat, r)
496    
497       while True:
498           normr=norm(r)
499           if normr<0: raise NegativeNorm
500           if verbose: print "iter: %s: norm(r) = %e, tolerance*norm(b) = %e"%(iter, normr,tolerance * normb)
501           if normr <= tolerance * normb: return x,r
502    
503           iter+=1 # k = iter = 1 first time through
504           if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
505    
506           q=Aprod(d)
507           alpha = rhat_dot_r / bilinearform(d, q)
508           x += alpha * d
509           r += (-alpha) * q
510    
511           rhat=Msolve(r)
512           rhat_dot_r_new = bilinearform(rhat, r)
513           beta = rhat_dot_r_new / rhat_dot_r
514           rhat+=beta * d
515           d=rhat
516    
517           rhat_dot_r = rhat_dot_r_new
518    
519    
520  class SaddlePointProblem(object):  class SaddlePointProblem(object):
521     """     """
522     This implements a solver for a saddlepoint problem     This implements a solver for a saddlepoint problem
# Line 621  class SaddlePointProblem(object): Line 746  class SaddlePointProblem(object):
746  #  #
747  #      return u,p  #      return u,p
748                        
749  # vim: expandtab shiftwidth=4:  def MaskFromBoundaryTag(function_space,*tags):
750       """
751       create a mask on the given function space which one for samples
752       that touch the boundary tagged by tags.
753    
754       usage: m=MaskFromBoundaryTag(Solution(domain),"left", "right")
755    
756       @param function_space: a given function space
757       @type function_space: L{escript.FunctionSpace}
758       @param tags: boundray tags
759       @type tags: C{str}
760       @return: a mask which marks samples used by C{function_space} that are touching the
761                boundary tagged by any of the given tags.
762       @rtype: L{escript.Data} of rank 0
763       """
764       pde=linearPDEs.LinearPDE(function_space.getDomain(),numEquations=1, numSolutions=1)
765       d=escript.Scalar(0.,escript.FunctionOnBoundary(function_space.getDomain()))
766       for t in tags: d.setTaggedValue(t,1.)
767       pde.setValue(y=d)
768       out=util.whereNonZero(pde.getRightHandSide())
769       if out.getFunctionSpace() == function_space:
770          return out
771       else:
772          return util.whereNonZero(util.interpolate(out,function_space))
773    

Legend:
Removed from v.1137  
changed lines
  Added in v.1312

  ViewVC Help
Powered by ViewVC 1.1.26