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

revision 1413 by trankine, Fri Jan 11 07:45:58 2008 UTC revision 1414 by gross, Thu Feb 14 10:01:43 2008 UTC
# Line 621  class ArithmeticTuple(object): Line 621  class ArithmeticTuple(object):
621             self.__items[i]+=other[i]             self.__items[i]+=other[i]
622         return self         return self
623
625          """
626          This provides a framwork for solving homogeneous saddle point problem of the form
627
628                 Av+B^*p=f
629                 Bv    =0
630
631          for the unknowns v and p and given operators A and B and given right hand side f.
632          B^* is the adjoint operator of B is the given inner product.
633
634          """
635          def __init__(self,**kwargs):
636            self.setTolerance()
637            self.setToleranceReductionFactor()
638
639          def initialize(self):
640            """
641            initialize the problem (overwrite)
642            """
643            pass
644          def B(self,v):
645             """
646             returns Bv (overwrite)
647             @rtype: equal to the type of p
648
649             @note: boundary conditions on p should be zero!
650             """
651             pass
652
653          def inner(self,p0,p1):
654             """
655             returns inner product of two element p0 and p1  (overwrite)
656
657             @type p0: equal to the type of p
658             @type p1: equal to the type of p
659             @rtype: C{float}
660
661             @rtype: equal to the type of p
662             """
663             pass
664
665          def solve_A(self,u,p):
666             """
667             solves Av=f-Au-B^*p with accuracy self.getReducedTolerance() (overwrite)
668
669             @rtype: equal to the type of v
670             @note: boundary conditions on v should be zero!
671             """
672             pass
673
674          def solve_prec(self,p):
675             """
676             provides a preconditioner for BA^{-1}B^* with accuracy self.getReducedTolerance() (overwrite)
677
678             @rtype: equal to the type of p
679             """
680             pass
681
682          def stoppingcriterium(self,Bv,v,p):
683             """
684             returns a True if iteration is terminated. (overwrite)
685
686             @rtype: C{bool}
687             """
688             pass
689
690          def __inner(self,p,r):
691             return self.inner(p,r[1])
692
693          def __stoppingcriterium(self,norm_r,r,p):
694              return self.stoppingcriterium(r[1],r[0],p)
695
696          def setTolerance(self,tolerance=1.e-8):
697                  self.__tol=tolerance
698          def getTolerance(self):
699                  return self.__tol
700          def setToleranceReductionFactor(self,reduction=0.01):
701                  self.__reduction=reduction
702          def getSubProblemTolerance(self):
703                  return self.__reduction*self.getTolerance()
704
705          def solve(self,v,p,max_iter=20, verbose=False, show_details=False):
706                  """
707                  solves the saddle point problem using initial guesses v and p.
708
709                  @param max_iter: maximum number of iteration steps.
710                  """
711                  self.verbose=verbose
712                  self.show_details=show_details and self.verbose
713
714              # Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask)
715
716              self.__z=v+self.solve_A(v,p*0)
717                  Bz=self.B(self.__z)
718                  #
719              #   solve BA^-1B^*p = Bz
720                  #
721                  #   note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv
722                  #
723                  #   with                    Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask)
724                  #                           A(v-z)=Az-B^*p-Az = f -Az - B^*p (v-z=0 on fixed_u_mask)
725                  #
726                  self.iter=0
727                  if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter
728                  p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p)
729              return r[0],p
730
731          def __Msolve(self,r):
732              return self.solve_prec(r[1])
733
734          def __Aprod(self,p):
735              # return BA^-1B*p
736              #solve Av =-B^*p as Av =f-Az-B^*p
737              v=self.solve_A(self.__z,p)
738              return ArithmeticTuple(v, self.B(v))
739
740