/[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 2264 by gross, Wed Feb 11 06:48:28 2009 UTC revision 2548 by jfenwick, Mon Jul 20 06:20:06 2009 UTC
# Line 1  Line 1 
1    
2  ########################################################  ########################################################
3  #  #
4  # Copyright (c) 2003-2008 by University of Queensland  # Copyright (c) 2003-2009 by University of Queensland
5  # Earth Systems Science Computational Center (ESSCC)  # Earth Systems Science Computational Center (ESSCC)
6  # http://www.uq.edu.au/esscc  # http://www.uq.edu.au/esscc
7  #  #
# Line 17  http://www.uq.edu.au/esscc Line 17  http://www.uq.edu.au/esscc
17  Primary Business: Queensland, Australia"""  Primary Business: Queensland, Australia"""
18  __license__="""Licensed under the Open Software License version 3.0  __license__="""Licensed under the Open Software License version 3.0
19  http://www.opensource.org/licenses/osl-3.0.php"""  http://www.opensource.org/licenses/osl-3.0.php"""
20  __url__="http://www.uq.edu.au/esscc/escript-finley"  __url__="https://launchpad.net/escript-finley"
21    
22  """  """
23  Provides some tools related to PDEs.  Provides some tools related to PDEs.
# Line 41  __author__="Lutz Gross, l.gross@uq.edu.a Line 41  __author__="Lutz Gross, l.gross@uq.edu.a
41    
42  import escript  import escript
43  import linearPDEs  import linearPDEs
44  import numarray  import numpy
45  import util  import util
46  import math  import math
47    
# Line 145  class Projector: Line 145  class Projector:
145      """      """
146      self.__pde = linearPDEs.LinearPDE(domain)      self.__pde = linearPDEs.LinearPDE(domain)
147      if fast:      if fast:
148          self.__pde.setSolverMethod(linearPDEs.LinearPDE.LUMPING)          self.__pde.getSolverOptions().setSolverMethod(linearPDEs.SolverOptions.LUMPING)
149      self.__pde.setSymmetryOn()      self.__pde.setSymmetryOn()
150      self.__pde.setReducedOrderTo(reduce)      self.__pde.setReducedOrderTo(reduce)
151      self.__pde.setValue(D = 1.)      self.__pde.setValue(D = 1.)
152      return      return
153      def getSolverOptions(self):
154        """
155        Returns the solver options of the PDE solver.
156        
157        @rtype: L{linearPDEs.SolverOptions}
158        """
159    
160    def __call__(self, input_data):    def __call__(self, input_data):
161      """      """
# Line 223  class NoPDE: Line 229  class NoPDE:
229           @param domain: domain of the PDE           @param domain: domain of the PDE
230           @type domain: L{Domain}           @type domain: L{Domain}
231           @param D: coefficient of the solution           @param D: coefficient of the solution
232           @type D: C{float}, C{int}, C{numarray.NumArray}, L{Data}           @type D: C{float}, C{int}, C{numpy.ndarray}, L{Data}
233           @param Y: right hand side           @param Y: right hand side
234           @type Y: C{float}, C{int}, C{numarray.NumArray}, L{Data}           @type Y: C{float}, C{int}, C{numpy.ndarray}, L{Data}
235           @param q: location of constraints           @param q: location of constraints
236           @type q: C{float}, C{int}, C{numarray.NumArray}, L{Data}           @type q: C{float}, C{int}, C{numpy.ndarray}, L{Data}
237           @param r: value of solution at locations of constraints           @param r: value of solution at locations of constraints
238           @type r: C{float}, C{int}, C{numarray.NumArray}, L{Data}           @type r: C{float}, C{int}, C{numpy.ndarray}, L{Data}
239           """           """
240           self.__domain=domain           self.__domain=domain
241           self.__D=D           self.__D=D
# Line 258  class NoPDE: Line 264  class NoPDE:
264           Assigns values to the parameters.           Assigns values to the parameters.
265    
266           @param D: coefficient of the solution           @param D: coefficient of the solution
267           @type D: C{float}, C{int}, C{numarray.NumArray}, L{Data}           @type D: C{float}, C{int}, C{numpy.ndarray}, L{Data}
268           @param Y: right hand side           @param Y: right hand side
269           @type Y: C{float}, C{int}, C{numarray.NumArray}, L{Data}           @type Y: C{float}, C{int}, C{numpy.ndarray}, L{Data}
270           @param q: location of constraints           @param q: location of constraints
271           @type q: C{float}, C{int}, C{numarray.NumArray}, L{Data}           @type q: C{float}, C{int}, C{numpy.ndarray}, L{Data}
272           @param r: value of solution at locations of constraints           @param r: value of solution at locations of constraints
273           @type r: C{float}, C{int}, C{numarray.NumArray}, L{Data}           @type r: C{float}, C{int}, C{numpy.ndarray}, L{Data}
274           """           """
275           if not D==None:           if not D==None:
276              self.__D=D              self.__D=D
# Line 312  class Locator: Line 318  class Locator:
318       given function space or domain which is closest to the given point x.       given function space or domain which is closest to the given point x.
319       """       """
320    
321       def __init__(self,where,x=numarray.zeros((3,))):       def __init__(self,where,x=numpy.zeros((3,))):
322         """         """
323         Initializes a Locator to access values in Data objects on the Doamin         Initializes a Locator to access values in Data objects on the Doamin
324         or FunctionSpace for the sample point which is closest to the given         or FunctionSpace for the sample point which is closest to the given
# Line 320  class Locator: Line 326  class Locator:
326    
327         @param where: function space         @param where: function space
328         @type where: L{escript.FunctionSpace}         @type where: L{escript.FunctionSpace}
329         @param x: coefficient of the solution         @param x: location(s) of the Locator
330         @type x: C{numarray.NumArray} or C{list} of C{numarray.NumArray}         @type x: C{numpy.ndarray} or C{list} of C{numpy.ndarray}
331         """         """
332         if isinstance(where,escript.FunctionSpace):         if isinstance(where,escript.FunctionSpace):
333            self.__function_space=where            self.__function_space=where
334         else:         else:
335            self.__function_space=escript.ContinuousFunction(where)            self.__function_space=escript.ContinuousFunction(where)
336           iterative=False
337         if isinstance(x, list):         if isinstance(x, list):
338               if len(x)==0:
339                  raise "ValueError", "At least one point must be given."
340               try:
341                 iter(x[0])
342                 iterative=True
343               except TypeError:
344                 iterative=False
345           if iterative:
346             self.__id=[]             self.__id=[]
347             for p in x:             for p in x:
348                self.__id.append(util.length(self.__function_space.getX()-p[:self.__function_space.getDim()]).minGlobalDataPoint())                self.__id.append(util.length(self.__function_space.getX()-p[:self.__function_space.getDim()]).minGlobalDataPoint())
# Line 390  class Locator: Line 405  class Locator:
405          object otherwise the object is returned.          object otherwise the object is returned.
406          """          """
407          if isinstance(data,escript.Data):          if isinstance(data,escript.Data):
408             if data.getFunctionSpace()==self.getFunctionSpace():             dat=util.interpolate(data,self.getFunctionSpace())
              dat=data  
            else:  
              dat=data.interpolate(self.getFunctionSpace())  
409             id=self.getId()             id=self.getId()
410             r=data.getRank()             r=data.getRank()
411             if isinstance(id,list):             if isinstance(id,list):
412                 out=[]                 out=[]
413                 for i in id:                 for i in id:
414                    o=data.getValueOfGlobalDataPoint(*i)                    o=numpy.array(dat.getTupleForGlobalDataPoint(*i))
415                    if data.getRank()==0:                    if data.getRank()==0:
416                       out.append(o[0])                       out.append(o[0])
417                    else:                    else:
418                       out.append(o)                       out.append(o)
419                 return out                 return out
420             else:             else:
421               out=data.getValueOfGlobalDataPoint(*id)               out=numpy.array(dat.getTupleForGlobalDataPoint(*id))
422               if data.getRank()==0:               if data.getRank()==0:
423                  return out[0]                  return out[0]
424               else:               else:
# Line 522  def PCG(r, Aprod, x, Msolve, bilinearfor Line 534  def PCG(r, Aprod, x, Msolve, bilinearfor
534         q=Aprod(d)         q=Aprod(d)
535         alpha = rhat_dot_r / bilinearform(d, q)         alpha = rhat_dot_r / bilinearform(d, q)
536         x += alpha * d         x += alpha * d
537         r += (-alpha) * q         if isinstance(q,ArithmeticTuple):
538           r += q * (-alpha)      # Doing it the other way calls the float64.__mul__ not AT.__rmul__
539           else:
540               r += (-alpha) * q
541         rhat=Msolve(r)         rhat=Msolve(r)
542         rhat_dot_r_new = bilinearform(rhat, r)         rhat_dot_r_new = bilinearform(rhat, r)
543         beta = rhat_dot_r_new / rhat_dot_r         beta = rhat_dot_r_new / rhat_dot_r
# Line 645  def NewtonGMRES(defect, x, iter_max=100, Line 659  def NewtonGMRES(defect, x, iter_max=100,
659     @type defect: L{Defect}     @type defect: L{Defect}
660     @param x: initial guess for the solution, C{x} is altered.     @param x: initial guess for the solution, C{x} is altered.
661     @type x: any object type allowing basic operations such as     @type x: any object type allowing basic operations such as
662              C{numarray.NumArray}, L{Data}              C{numpy.ndarray}, L{Data}
663     @param iter_max: maximum number of iteration steps     @param iter_max: maximum number of iteration steps
664     @type iter_max: positive C{int}     @type iter_max: positive C{int}
665     @param sub_iter_max: maximum number of inner iteration steps     @param sub_iter_max: maximum number of inner iteration steps
# Line 728  def __givapp(c,s,vin): Line 742  def __givapp(c,s,vin):
742          for i in range(len(c)):          for i in range(len(c)):
743              w1=c[i]*vrot[i]-s[i]*vrot[i+1]              w1=c[i]*vrot[i]-s[i]*vrot[i+1]
744          w2=s[i]*vrot[i]+c[i]*vrot[i+1]          w2=s[i]*vrot[i]+c[i]*vrot[i+1]
745              vrot[i:i+2]=w1,w2              vrot[i]=w1
746                vrot[i+1]=w2
747      return vrot      return vrot
748    
749  def __FDGMRES(F0, defect, x0, atol, iter_max=100, iter_restart=20):  def __FDGMRES(F0, defect, x0, atol, iter_max=100, iter_restart=20):
750     h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)     h=numpy.zeros((iter_restart,iter_restart),numpy.float64)
751     c=numarray.zeros(iter_restart,numarray.Float64)     c=numpy.zeros(iter_restart,numpy.float64)
752     s=numarray.zeros(iter_restart,numarray.Float64)     s=numpy.zeros(iter_restart,numpy.float64)
753     g=numarray.zeros(iter_restart,numarray.Float64)     g=numpy.zeros(iter_restart,numpy.float64)
754     v=[]     v=[]
755    
756     rho=defect.norm(F0)     rho=defect.norm(F0)
# Line 776  def __FDGMRES(F0, defect, x0, atol, iter Line 791  def __FDGMRES(F0, defect, x0, atol, iter
791    
792          #   Form and store the information for the new Givens rotation          #   Form and store the information for the new Givens rotation
793          if iter > 0 :          if iter > 0 :
794              hhat=numarray.zeros(iter+1,numarray.Float64)              hhat=numpy.zeros(iter+1,numpy.float64)
795              for i in range(iter+1) : hhat[i]=h[i,iter]              for i in range(iter+1) : hhat[i]=h[i,iter]
796              hhat=__givapp(c[0:iter],s[0:iter],hhat);              hhat=__givapp(c[0:iter],s[0:iter],hhat);
797              for i in range(iter+1) : h[i,iter]=hhat[i]              for i in range(iter+1) : h[i,iter]=hhat[i]
# Line 788  def __FDGMRES(F0, defect, x0, atol, iter Line 803  def __FDGMRES(F0, defect, x0, atol, iter
803              s[iter]=-h[iter+1,iter]/mu              s[iter]=-h[iter+1,iter]/mu
804              h[iter,iter]=c[iter]*h[iter,iter]-s[iter]*h[iter+1,iter]              h[iter,iter]=c[iter]*h[iter,iter]-s[iter]*h[iter+1,iter]
805              h[iter+1,iter]=0.0              h[iter+1,iter]=0.0
806              g[iter:iter+2]=__givapp(c[iter],s[iter],g[iter:iter+2])              gg=__givapp(c[iter],s[iter],[g[iter],g[iter+1]])
807                g[iter]=gg[0]
808                g[iter+1]=gg[1]
809    
810          # Update the residual norm          # Update the residual norm
811          rho=abs(g[iter+1])          rho=abs(g[iter+1])
# Line 797  def __FDGMRES(F0, defect, x0, atol, iter Line 814  def __FDGMRES(F0, defect, x0, atol, iter
814     # At this point either iter > iter_max or rho < tol.     # At this point either iter > iter_max or rho < tol.
815     # It's time to compute x and leave.     # It's time to compute x and leave.
816     if iter > 0 :     if iter > 0 :
817       y=numarray.zeros(iter,numarray.Float64)       y=numpy.zeros(iter,numpy.float64)
818       y[iter-1] = g[iter-1] / h[iter-1,iter-1]       y[iter-1] = g[iter-1] / h[iter-1,iter-1]
819       if iter > 1 :       if iter > 1 :
820          i=iter-2          i=iter-2
821          while i>=0 :          while i>=0 :
822            y[i] = ( g[i] - numarray.dot(h[i,i+1:iter], y[i+1:iter])) / h[i,i]            y[i] = ( g[i] - numpy.dot(h[i,i+1:iter], y[i+1:iter])) / h[i,i]
823            i=i-1            i=i-1
824       xhat=v[iter-1]*y[iter-1]       xhat=v[iter-1]*y[iter-1]
825       for i in range(iter-1):       for i in range(iter-1):
# Line 892  def GMRES(r, Aprod, x, bilinearform, ato Line 909  def GMRES(r, Aprod, x, bilinearform, ato
909  def _GMRESm(r, Aprod, x, bilinearform, atol, iter_max=100, iter_restart=20, verbose=False, P_R=None):  def _GMRESm(r, Aprod, x, bilinearform, atol, iter_max=100, iter_restart=20, verbose=False, P_R=None):
910     iter=0     iter=0
911    
912     h=numarray.zeros((iter_restart+1,iter_restart),numarray.Float64)     h=numpy.zeros((iter_restart+1,iter_restart),numpy.float64)
913     c=numarray.zeros(iter_restart,numarray.Float64)     c=numpy.zeros(iter_restart,numpy.float64)
914     s=numarray.zeros(iter_restart,numarray.Float64)     s=numpy.zeros(iter_restart,numpy.float64)
915     g=numarray.zeros(iter_restart+1,numarray.Float64)     g=numpy.zeros(iter_restart+1,numpy.float64)
916     v=[]     v=[]
917    
918     r_dot_r = bilinearform(r, r)     r_dot_r = bilinearform(r, r)
# Line 949  def _GMRESm(r, Aprod, x, bilinearform, a Line 966  def _GMRESm(r, Aprod, x, bilinearform, a
966          s[iter]=-h[iter+1,iter]/mu          s[iter]=-h[iter+1,iter]/mu
967          h[iter,iter]=c[iter]*h[iter,iter]-s[iter]*h[iter+1,iter]          h[iter,iter]=c[iter]*h[iter,iter]-s[iter]*h[iter+1,iter]
968          h[iter+1,iter]=0.0          h[iter+1,iter]=0.0
969          g[iter:iter+2]=__givapp(c[iter],s[iter],g[iter:iter+2])                  gg=__givapp(c[iter],s[iter],[g[iter],g[iter+1]])
970                    g[iter]=gg[0]
971                    g[iter+1]=gg[1]
972  # Update the residual norm  # Update the residual norm
973    
974          rho=abs(g[iter+1])          rho=abs(g[iter+1])
# Line 961  def _GMRESm(r, Aprod, x, bilinearform, a Line 980  def _GMRESm(r, Aprod, x, bilinearform, a
980    
981     if verbose: print "GMRES: iteration stopped after %s step."%iter     if verbose: print "GMRES: iteration stopped after %s step."%iter
982     if iter > 0 :     if iter > 0 :
983       y=numarray.zeros(iter,numarray.Float64)       y=numpy.zeros(iter,numpy.float64)
984       y[iter-1] = g[iter-1] / h[iter-1,iter-1]       y[iter-1] = g[iter-1] / h[iter-1,iter-1]
985       if iter > 1 :       if iter > 1 :
986          i=iter-2          i=iter-2
987          while i>=0 :          while i>=0 :
988            y[i] = ( g[i] - numarray.dot(h[i,i+1:iter], y[i+1:iter])) / h[i,i]            y[i] = ( g[i] - numpy.dot(h[i,i+1:iter], y[i+1:iter])) / h[i,i]
989            i=i-1            i=i-1
990       xhat=v[iter-1]*y[iter-1]       xhat=v[iter-1]*y[iter-1]
991       for i in range(iter-1):       for i in range(iter-1):
# Line 1274  class ArithmeticTuple(object): Line 1293  class ArithmeticTuple(object):
1293     Example of usage::     Example of usage::
1294    
1295         from esys.escript import Data         from esys.escript import Data
1296         from numarray import array         from numpy import array
1297         a=Data(...)         a=Data(...)
1298         b=array([1.,4.])         b=array([1.,4.])
1299         x=ArithmeticTuple(a,b)         x=ArithmeticTuple(a,b)
# Line 1458  class HomogeneousSaddlePointProblem(obje Line 1477  class HomogeneousSaddlePointProblem(obje
1477        for the unknowns M{v} and M{p} and given operators M{A} and M{B} and        for the unknowns M{v} and M{p} and given operators M{A} and M{B} and
1478        given right hand side M{f}. M{B^*} is the adjoint operator of M{B}.        given right hand side M{f}. M{B^*} is the adjoint operator of M{B}.
1479        """        """
1480        def __init__(self,**kwargs):        def __init__(self, adaptSubTolerance=True, **kwargs):
1481        """
1482        initializes the saddle point problem
1483        
1484        @param adaptSubTolerance: If True the tolerance for subproblem is set automatically.
1485        @type adaptSubTolerance: C{bool}
1486        """
1487          self.setTolerance()          self.setTolerance()
1488          self.setAbsoluteTolerance()          self.setAbsoluteTolerance()
1489          self.setSubProblemTolerance()      self.__adaptSubTolerance=adaptSubTolerance
   
1490        #=============================================================        #=============================================================
1491        def initialize(self):        def initialize(self):
1492          """          """
# Line 1470  class HomogeneousSaddlePointProblem(obje Line 1494  class HomogeneousSaddlePointProblem(obje
1494          """          """
1495          pass          pass
1496    
1497        def inner_pBv(self,p,v):        def inner_pBv(self,p,Bv):
1498           """           """
1499           Returns inner product of element p and Bv (overwrite).           Returns inner product of element p and Bv (overwrite).
1500    
# Line 1480  class HomogeneousSaddlePointProblem(obje Line 1504  class HomogeneousSaddlePointProblem(obje
1504           @rtype: C{float}           @rtype: C{float}
1505           @note: used if PCG is applied.           @note: used if PCG is applied.
1506           """           """
1507           raise NotImplementedError,"no inner product for p implemented."           raise NotImplementedError,"no inner product for p and Bv implemented."
1508    
1509        def inner_p(self,p0,p1):        def inner_p(self,p0,p1):
1510           """           """
# Line 1502  class HomogeneousSaddlePointProblem(obje Line 1526  class HomogeneousSaddlePointProblem(obje
1526           @rtype: non-negative C{float}           @rtype: non-negative C{float}
1527           """           """
1528           raise NotImplementedError,"no norm of v implemented."           raise NotImplementedError,"no norm of v implemented."
   
   
1529        def getV(self, p, v0):        def getV(self, p, v0):
1530           """           """
1531           return the value for v for a given p (overwrite)           return the value for v for a given p (overwrite)
# Line 1513  class HomogeneousSaddlePointProblem(obje Line 1535  class HomogeneousSaddlePointProblem(obje
1535           @return: v given as M{v= A^{-1} (f-B^*p)}           @return: v given as M{v= A^{-1} (f-B^*p)}
1536           """           """
1537           raise NotImplementedError,"no v calculation implemented."           raise NotImplementedError,"no v calculation implemented."
1538    
1539                    
1540        def norm_Bv(self,v):        def Bv(self,v):
1541          """          """
1542          Returns Bv (overwrite).          Returns Bv (overwrite).
1543    
# Line 1523  class HomogeneousSaddlePointProblem(obje Line 1546  class HomogeneousSaddlePointProblem(obje
1546          """          """
1547          raise NotImplementedError, "no operator B implemented."          raise NotImplementedError, "no operator B implemented."
1548    
1549          def norm_Bv(self,Bv):
1550            """
1551            Returns the norm of Bv (overwrite).
1552    
1553            @rtype: equal to the type of p
1554            @note: boundary conditions on p should be zero!
1555            """
1556            raise NotImplementedError, "no norm of Bv implemented."
1557    
1558        def solve_AinvBt(self,p):        def solve_AinvBt(self,p):
1559           """           """
1560           Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}           Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}
# Line 1534  class HomogeneousSaddlePointProblem(obje Line 1566  class HomogeneousSaddlePointProblem(obje
1566           """           """
1567           raise NotImplementedError,"no operator A implemented."           raise NotImplementedError,"no operator A implemented."
1568    
1569        def solve_precB(self,v):        def solve_prec(self,Bv):
1570           """           """
1571           Provides a preconditioner for M{BA^{-1}B^*} with accuracy           Provides a preconditioner for M{BA^{-1}B^*} applied to Bv with accuracy
1572           L{self.getSubProblemTolerance()} (overwrite).           L{self.getSubProblemTolerance()} (overwrite).
1573    
1574           @rtype: equal to the type of p           @rtype: equal to the type of p
1575           @note: boundary conditions on p should be zero!           @note: boundary conditions on p should be zero!
1576           """           """
1577           raise NotImplementedError,"no preconditioner for Schur complement implemented."           raise NotImplementedError,"no preconditioner for Schur complement implemented."
1578          def setSubProblemTolerance(self):
1579             """
1580         Updates the tolerance for subproblems
1581         @note: method is typically the method is overwritten.
1582             """
1583             pass
1584        #=============================================================        #=============================================================
1585        def __Aprod_PCG(self,p):        def __Aprod_PCG(self,p):
1586            return self.solve_AinvBt(p)            dv=self.solve_AinvBt(p)
1587              return ArithmeticTuple(dv,self.Bv(dv))
1588    
1589        def __inner_PCG(self,p,v):        def __inner_PCG(self,p,r):
1590           return self.inner_pBv(p,v)           return self.inner_pBv(p,r[1])
1591    
1592        def __Msolve_PCG(self,v):        def __Msolve_PCG(self,r):
1593            return self.solve_precB(v)            return self.solve_prec(r[1])
1594        #=============================================================        #=============================================================
1595        def __Aprod_GMRES(self,p):        def __Aprod_GMRES(self,p):
1596            return self.solve_precB(self.solve_AinvBt(p))            return self.solve_prec(self.Bv(self.solve_AinvBt(p)))
1597        def __inner_GMRES(self,p0,p1):        def __inner_GMRES(self,p0,p1):
1598           return self.inner_p(p0,p1)           return self.inner_p(p0,p1)
1599    
1600        #=============================================================        #=============================================================
1601        def norm_p(self,p):        def norm_p(self,p):
1602            """            """
# Line 1569  class HomogeneousSaddlePointProblem(obje Line 1609  class HomogeneousSaddlePointProblem(obje
1609            f=self.inner_p(p,p)            f=self.inner_p(p,p)
1610            if f<0: raise ValueError,"negative pressure norm."            if f<0: raise ValueError,"negative pressure norm."
1611            return math.sqrt(f)            return math.sqrt(f)
1612                    def adaptSubTolerance(self):
1613          """
1614        def solve(self,v,p,max_iter=20, verbose=False, show_details=False, usePCG=True, iter_restart=20, max_correction_steps=10):        Returns True if tolerance adaption for subproblem is choosen.
1615          """
1616              self.__adaptSubTolerance
1617          
1618          def solve(self,v,p,max_iter=20, verbose=False, usePCG=True, iter_restart=20, max_correction_steps=10):
1619           """           """
1620           Solves the saddle point problem using initial guesses v and p.           Solves the saddle point problem using initial guesses v and p.
1621    
# Line 1584  class HomogeneousSaddlePointProblem(obje Line 1628  class HomogeneousSaddlePointProblem(obje
1628                            attempt                            attempt
1629           @param verbose: if True, shows information on the progress of the           @param verbose: if True, shows information on the progress of the
1630                           saddlepoint problem solver.                           saddlepoint problem solver.
          @param show_details: if True, shows details of the sub problem solver  
1631           @param iter_restart: restart the iteration after C{iter_restart} steps           @param iter_restart: restart the iteration after C{iter_restart} steps
1632                                (only used if useUzaw=False)                                (only used if useUzaw=False)
1633           @type usePCG: C{bool}           @type usePCG: C{bool}
1634           @type max_iter: C{int}           @type max_iter: C{int}
1635           @type verbose: C{bool}           @type verbose: C{bool}
          @type show_details: C{bool}  
1636           @type iter_restart: C{int}           @type iter_restart: C{int}
1637           @rtype: C{tuple} of L{Data} objects           @rtype: C{tuple} of L{Data} objects
1638           """           """
1639           self.verbose=verbose           self.verbose=verbose
          self.show_details=show_details and self.verbose  
1640           rtol=self.getTolerance()           rtol=self.getTolerance()
1641           atol=self.getAbsoluteTolerance()           atol=self.getAbsoluteTolerance()
1642         if self.adaptSubTolerance(): self.setSubProblemTolerance()
1643           correction_step=0           correction_step=0
1644           converged=False           converged=False
1645           while not converged:           while not converged:
1646                # calculate velocity for current pressure:                # calculate velocity for current pressure:
1647                v=self.getV(p,v)                v=self.getV(p,v)
1648                #                Bv=self.Bv(v)
1649                norm_v=self.norm_v(v)                norm_v=self.norm_v(v)
1650                norm_Bv=self.norm_Bv(v)                norm_Bv=self.norm_Bv(Bv)
1651                ATOL=norm_v*rtol+atol                ATOL=norm_v*rtol+atol
1652                if self.verbose: print "saddle point solver: norm v= %e, norm_Bv= %e, tolerance = %e."%(norm_v, norm_Bv,ATOL)                if self.verbose: print "HomogeneousSaddlePointProblem: norm v= %e, norm_Bv= %e, tolerance = %e."%(norm_v, norm_Bv,ATOL)
1653                if not ATOL>0: raise ValueError,"overall absolute tolerance needs to be positive."                if not ATOL>0: raise ValueError,"overall absolute tolerance needs to be positive."
1654                if norm_Bv <= ATOL:                if norm_Bv <= ATOL:
1655                   converged=True                   converged=True
# Line 1615  class HomogeneousSaddlePointProblem(obje Line 1657  class HomogeneousSaddlePointProblem(obje
1657                   correction_step+=1                   correction_step+=1
1658                   if correction_step>max_correction_steps:                   if correction_step>max_correction_steps:
1659                        raise CorrectionFailed,"Given up after %d correction steps."%correction_step                        raise CorrectionFailed,"Given up after %d correction steps."%correction_step
1660                   dp=self.solve_precB(v)                   dp=self.solve_prec(Bv)
1661                   if usePCG:                   if usePCG:
1662                     norm2=self.inner_pBv(dp,v)                     norm2=self.inner_pBv(dp,Bv)
1663                     if norm2<0: raise ValueError,"negative PCG norm."                     if norm2<0: raise ValueError,"negative PCG norm."
1664                     norm2=math.sqrt(norm2)                     norm2=math.sqrt(norm2)
1665                   else:                   else:
1666                     norm2=self.norm_p(dp)                     norm2=self.norm_p(dp)
1667                   ATOL_ITER=ATOL/norm_Bv*norm2                   ATOL_ITER=ATOL/norm_Bv*norm2*0.5
1668                   if self.verbose: print "saddle point solver: tolerance for solver: %e"%ATOL_ITER                   if self.verbose: print "HomogeneousSaddlePointProblem: tolerance for solver: %e"%ATOL_ITER
1669                   if usePCG:                   if usePCG:
1670                         p,v0,a_norm=PCG(v,self.__Aprod_PCG,p,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL_ITER, rtol=0.,iter_max=max_iter, verbose=self.verbose)                         p,v0,a_norm=PCG(ArithmeticTuple(v,Bv),self.__Aprod_PCG,p,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL_ITER, rtol=0.,iter_max=max_iter, verbose=self.verbose)
1671                   else:                   else:
1672                         p=GMRES(dp,self.__Aprod_GMRES, p, self.__inner_GMRES,atol=ATOL_ITER, rtol=0.,iter_max=max_iter, iter_restart=iter_restart, verbose=self.verbose)                         p=GMRES(dp,self.__Aprod_GMRES, p, self.__inner_GMRES,atol=ATOL_ITER, rtol=0.,iter_max=max_iter, iter_restart=iter_restart, verbose=self.verbose)
1673           if self.verbose: print "saddle point solver: tolerance reached."           if self.verbose: print "HomogeneousSaddlePointProblem: tolerance reached."
1674       return v,p       return v,p
1675    
1676        #========================================================================        #========================================================================
# Line 1642  class HomogeneousSaddlePointProblem(obje Line 1684  class HomogeneousSaddlePointProblem(obje
1684           if tolerance<0:           if tolerance<0:
1685               raise ValueError,"tolerance must be positive."               raise ValueError,"tolerance must be positive."
1686           self.__rtol=tolerance           self.__rtol=tolerance
          self.setSubProblemTolerance()  
1687    
1688        def getTolerance(self):        def getTolerance(self):
1689           """           """
# Line 1673  class HomogeneousSaddlePointProblem(obje Line 1714  class HomogeneousSaddlePointProblem(obje
1714           """           """
1715           return self.__atol           return self.__atol
1716    
1717        def setSubProblemTolerance(self,rtol=None):        def getSubProblemTolerance(self):
1718           """           """
1719           Sets the relative tolerance to solve the subproblem(s).           Sets the relative tolerance to solve the subproblem(s).
1720    
1721           @param rtol: relative tolerence           @param rtol: relative tolerence
1722           @type rtol: positive C{float}           @type rtol: positive C{float}
1723           """           """
1724           if rtol == None:           return max(200.*util.EPSILON,self.getTolerance()**2)
               rtol=max(200.*util.EPSILON,self.getTolerance()**2)  
          if rtol<=0:  
              raise ValueError,"tolerance must be positive."  
          self.__sub_tol=rtol  
   
       def getSubProblemTolerance(self):  
          """  
          Returns the subproblem reduction factor.  
   
          @return: subproblem reduction factor  
          @rtype: C{float}  
          """  
          return self.__sub_tol  
1725    
1726  def MaskFromBoundaryTag(domain,*tags):  def MaskFromBoundaryTag(domain,*tags):
1727     """     """
# Line 1716  def MaskFromBoundaryTag(domain,*tags): Line 1744  def MaskFromBoundaryTag(domain,*tags):
1744     pde.setValue(y=d)     pde.setValue(y=d)
1745     return util.whereNonZero(pde.getRightHandSide())     return util.whereNonZero(pde.getRightHandSide())
1746    
1747  #==============================================================================  def MaskFromTag(domain,*tags):
 class SaddlePointProblem(object):  
1748     """     """
1749     This implements a solver for a saddle point problem     Creates a mask on the Solution(domain) function space where the value is
1750       one for samples that touch regions tagged by tags.
1751    
1752     M{f(u,p)=0}     Usage: m=MaskFromTag(domain, "ham")
1753     M{g(u)=0}  
1754       @param domain: domain to be used
1755       @type domain: L{escript.Domain}
1756       @param tags: boundary tags
1757       @type tags: C{str}
1758       @return: a mask which marks samples that are touching the boundary tagged
1759                by any of the given tags
1760       @rtype: L{escript.Data} of rank 0
1761       """
1762       pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1763       d=escript.Scalar(0.,escript.Function(domain))
1764       for t in tags: d.setTaggedValue(t,1.)
1765       pde.setValue(Y=d)
1766       return util.whereNonZero(pde.getRightHandSide())
1767    
    for u and p. The problem is solved with an inexact Uszawa scheme for p:  
   
    M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})}  
    M{Q_g (p^{k+1}-p^{k}) =   g(u^{k+1})}  
   
    where Q_f is an approximation of the Jacobian A_f of f with respect to u and  
    Q_f is an approximation of A_g A_f^{-1} A_g with A_g is the Jacobian of g  
    with respect to p. As a the construction of a 'proper' Q_g can be difficult,  
    non-linear conjugate gradient method is applied to solve for p, so Q_g plays  
    in fact the role of a preconditioner.  
    """  
    def __init__(self,verbose=False,*args):  
        """  
        Initializes the problem.  
   
        @param verbose: if True, some information is printed in the process  
        @type verbose: C{bool}  
        @note: this method may be overwritten by a particular saddle point  
               problem  
        """  
        print "SaddlePointProblem should not be used anymore!"  
        if not isinstance(verbose,bool):  
             raise TypeError("verbose needs to be of type bool.")  
        self.__verbose=verbose  
        self.relaxation=1.  
        DeprecationWarning("SaddlePointProblem should not be used anymore.")  
   
    def trace(self,text):  
        """  
        Prints C{text} if verbose has been set.  
   
        @param text: a text message  
        @type text: C{str}  
        """  
        if self.__verbose: print "%s: %s"%(str(self),text)  
   
    def solve_f(self,u,p,tol=1.e-8):  
        """  
        Solves  
   
        A_f du = f(u,p)  
   
        with tolerance C{tol} and returns du. A_f is Jacobian of f with respect  
        to u.  
   
        @param u: current approximation of u  
        @type u: L{escript.Data}  
        @param p: current approximation of p  
        @type p: L{escript.Data}  
        @param tol: tolerance expected for du  
        @type tol: C{float}  
        @return: increment du  
        @rtype: L{escript.Data}  
        @note: this method has to be overwritten by a particular saddle point  
               problem  
        """  
        pass  
   
    def solve_g(self,u,tol=1.e-8):  
        """  
        Solves  
   
        Q_g dp = g(u)  
   
        where Q_g is a preconditioner for A_g A_f^{-1} A_g with A_g is the  
        Jacobian of g with respect to p.  
   
        @param u: current approximation of u  
        @type u: L{escript.Data}  
        @param tol: tolerance expected for dp  
        @type tol: C{float}  
        @return: increment dp  
        @rtype: L{escript.Data}  
        @note: this method has to be overwritten by a particular saddle point  
               problem  
        """  
        pass  
   
    def inner(self,p0,p1):  
        """  
        Inner product of p0 and p1 approximating p. Typically this returns  
        C{integrate(p0*p1)}.  
        @return: inner product of p0 and p1  
        @rtype: C{float}  
        """  
        pass  
   
    subiter_max=3  
    def solve(self,u0,p0,tolerance=1.e-6,tolerance_u=None,iter_max=100,accepted_reduction=0.995,relaxation=None):  
         """  
         Runs the solver.  
   
         @param u0: initial guess for C{u}  
         @type u0: L{esys.escript.Data}  
         @param p0: initial guess for C{p}  
         @type p0: L{esys.escript.Data}  
         @param tolerance: tolerance for relative error in C{u} and C{p}  
         @type tolerance: positive C{float}  
         @param tolerance_u: tolerance for relative error in C{u} if different  
                             from C{tolerance}  
         @type tolerance_u: positive C{float}  
         @param iter_max: maximum number of iteration steps  
         @type iter_max: C{int}  
         @param accepted_reduction: if the norm g cannot be reduced by  
                                    C{accepted_reduction} backtracking to adapt  
                                    the relaxation factor. If  
                                    C{accepted_reduction=None} no backtracking  
                                    is used.  
         @type accepted_reduction: positive C{float} or C{None}  
         @param relaxation: initial relaxation factor. If C{relaxation==None},  
                            the last relaxation factor is used.  
         @type relaxation: C{float} or C{None}  
         """  
         tol=1.e-2  
         if tolerance_u==None: tolerance_u=tolerance  
         if not relaxation==None: self.relaxation=relaxation  
         if accepted_reduction ==None:  
               angle_limit=0.  
         elif accepted_reduction>=1.:  
               angle_limit=0.  
         else:  
               angle_limit=util.sqrt(1-accepted_reduction**2)  
         self.iter=0  
         u=u0  
         p=p0  
         #  
         #   initialize things:  
         #  
         converged=False  
         #  
         #  start loop:  
         #  
         #  initial search direction is g  
         #  
         while not converged :  
             if self.iter>iter_max:  
                 raise ArithmeticError("no convergence after %s steps."%self.iter)  
             f_new=self.solve_f(u,p,tol)  
             norm_f_new = util.Lsup(f_new)  
             u_new=u-f_new  
             g_new=self.solve_g(u_new,tol)  
             self.iter+=1  
             norm_g_new = util.sqrt(self.inner(g_new,g_new))  
             if norm_f_new==0. and norm_g_new==0.: return u, p  
             if self.iter>1 and not accepted_reduction==None:  
                #  
                #   did we manage to reduce the norm of G? I  
                #   if not we start a backtracking procedure  
                #  
                # print "new/old norm = ",norm_g_new, norm_g, norm_g_new/norm_g  
                if norm_g_new > accepted_reduction * norm_g:  
                   sub_iter=0  
                   s=self.relaxation  
                   d=g  
                   g_last=g  
                   self.trace("    start substepping: f = %s, g = %s, relaxation = %s."%(norm_f_new, norm_g_new, s))  
                   while sub_iter < self.subiter_max and  norm_g_new > accepted_reduction * norm_g:  
                      dg= g_new-g_last  
                      norm_dg=abs(util.sqrt(self.inner(dg,dg))/self.relaxation)  
                      rad=self.inner(g_new,dg)/self.relaxation  
                      # print "   ",sub_iter,": rad, norm_dg:",abs(rad), norm_dg*norm_g_new * angle_limit  
                      # print "   ",sub_iter,": rad, norm_dg:",rad, norm_dg, norm_g_new, norm_g  
                      if abs(rad) < norm_dg*norm_g_new * angle_limit:  
                          if sub_iter>0: self.trace("    no further improvements expected from backtracking.")  
                          break  
                      r=self.relaxation  
                      self.relaxation= - rad/norm_dg**2  
                      s+=self.relaxation  
                      #####  
                      # a=g_new+self.relaxation*dg/r  
                      # print "predicted new norm = ",util.sqrt(self.inner(a,a)),util.sqrt(self.inner(g_new,g_new)), self.relaxation  
                      #####  
                      g_last=g_new  
                      p+=self.relaxation*d  
                      f_new=self.solve_f(u,p,tol)  
                      u_new=u-f_new  
                      g_new=self.solve_g(u_new,tol)  
                      self.iter+=1  
                      norm_f_new = util.Lsup(f_new)  
                      norm_g_new = util.sqrt(self.inner(g_new,g_new))  
                      # print "   ",sub_iter," new g norm",norm_g_new  
                      self.trace("    %s th sub-step: f = %s, g = %s, relaxation = %s."%(sub_iter, norm_f_new, norm_g_new, s))  
                      #  
                      #   can we expect reduction of g?  
                      #  
                      # u_last=u_new  
                      sub_iter+=1  
                   self.relaxation=s  
             #  
             #  check for convergence:  
             #  
             norm_u_new = util.Lsup(u_new)  
             p_new=p+self.relaxation*g_new  
             norm_p_new = util.sqrt(self.inner(p_new,p_new))  
             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))  
   
             if self.iter>1:  
                dg2=g_new-g  
                df2=f_new-f  
                norm_dg2=util.sqrt(self.inner(dg2,dg2))  
                norm_df2=util.Lsup(df2)  
                # print norm_g_new, norm_g, norm_dg, norm_p, tolerance  
                tol_eq_g=tolerance*norm_dg2/(norm_g*abs(self.relaxation))*norm_p_new  
                tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new  
                if norm_g_new <= tol_eq_g and norm_f_new <= tol_eq_f:  
                    converged=True  
             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  
         self.trace("convergence after %s steps."%self.iter)  
         return u,p  
1768    

Legend:
Removed from v.2264  
changed lines
  Added in v.2548

  ViewVC Help
Powered by ViewVC 1.1.26