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

revision 2245 by gross, Wed Feb 4 06:27:59 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 534  def PCG(r, Aprod, x, Msolve, bilinearfor Line 548  def PCG(r, Aprod, x, Msolve, bilinearfor
548         if rhat_dot_r<0: raise NegativeNorm,"negative norm."         if rhat_dot_r<0: raise NegativeNorm,"negative norm."
549         if verbose: print "PCG: iteration step %s: residual norm = %e"%(iter, math.sqrt(rhat_dot_r))         if verbose: print "PCG: iteration step %s: residual norm = %e"%(iter, math.sqrt(rhat_dot_r))
550     if verbose: print "PCG: tolerance reached after %s steps."%iter     if verbose: print "PCG: tolerance reached after %s steps."%iter
551     return x,r     return x,r,math.sqrt(rhat_dot_r)
552
553  class Defect(object):  class Defect(object):
554      """      """
# 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 817  def __FDGMRES(F0, defect, x0, atol, iter Line 834  def __FDGMRES(F0, defect, x0, atol, iter
834
835     return xhat,stopped     return xhat,stopped
836
837  def GMRES(r, Aprod, x, bilinearform, atol=0, rtol=1.e-8, iter_max=100, iter_restart=20, verbose=False):  def GMRES(r, Aprod, x, bilinearform, atol=0, rtol=1.e-8, iter_max=100, iter_restart=20, verbose=False,P_R=None):
838     """     """
839     Solver for     Solver for
840
# Line 881  def GMRES(r, Aprod, x, bilinearform, ato Line 898  def GMRES(r, Aprod, x, bilinearform, ato
898        else:        else:
899           r2=1*r           r2=1*r
900        x2=x*1.        x2=x*1.
901        x,stopped=_GMRESm(r2, Aprod, x, bilinearform, atol2, iter_max=iter_max-iter, iter_restart=m, verbose=verbose)        x,stopped=_GMRESm(r2, Aprod, x, bilinearform, atol2, iter_max=iter_max-iter, iter_restart=m, verbose=verbose,P_R=P_R)
902        iter+=iter_restart        iter+=iter_restart
903        if stopped: break        if stopped: break
904        if verbose: print "GMRES: restart."        if verbose: print "GMRES: restart."
905        restarted=True        restarted=True
906     if verbose: print "GMRES: tolerance has reached."     if verbose: print "GMRES: tolerance has been reached."
907     return x     return x
908
909  def _GMRESm(r, Aprod, x, bilinearform, atol, iter_max=100, iter_restart=20, verbose=False):  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 910  def _GMRESm(r, Aprod, x, bilinearform, a Line 927  def _GMRESm(r, Aprod, x, bilinearform, a
927
928      if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max      if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
929
930      p=Aprod(v[iter])          if P_R!=None:
931                p=Aprod(P_R(v[iter]))
932            else:
933            p=Aprod(v[iter])
934      v.append(p)      v.append(p)
935
936      v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))      v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
# Line 946  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 958  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):
992      xhat += v[i]*y[i]      xhat += v[i]*y[i]
993     else:     else:
994       xhat=v[0] * 0       xhat=v[0] * 0
995       if P_R!=None:
996     x += xhat        x += P_R(xhat)
997       else:
998          x += xhat
999     if iter<iter_restart-1:     if iter<iter_restart-1:
1000        stopped=True        stopped=True
1001     else:     else:
# Line 1269  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 1453  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 1465  class HomogeneousSaddlePointProblem(obje Line 1494  class HomogeneousSaddlePointProblem(obje
1494          """          """
1495          pass          pass
1496
def B(self,v):
"""
Returns Bv (overwrite).

@rtype: equal to the type of p
@note: boundary conditions on p should be zero!
"""
raise NotImplementedError, "no operator B implemented."

1497        def inner_pBv(self,p,Bv):        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
1501           @type p: equal to the type of p           @param p: a pressure increment
1502           @type Bv: equal to the type of result of operator B           @param v: a residual
1503           @rtype: equal to the type of p           @return: inner product of element p and Bv
1504             @rtype: C{float}
1505             @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           """           """
1511           Returns inner product of element p0 and p1 (overwrite).           Returns inner product of p0 and p1 (overwrite).
1512
1513           @type p0: equal to the type of p           @param p0: a pressure
1514           @type p1: equal to the type of p           @param p1: a pressure
1515           @rtype: equal to the type of p           @return: inner product of p0 and p1
1516             @rtype: C{float}
1517           """           """
1518           raise NotImplementedError,"no inner product for p implemented."           raise NotImplementedError,"no inner product for p implemented."
1519
1520          def norm_v(self,v):
1521             """
1522             Returns the norm of v (overwrite).
1523
1524        def inner_v(self,v0,v1):           @param v: a velovity
1525             @return: norm of v
1526             @rtype: non-negative C{float}
1527             """
1528             raise NotImplementedError,"no norm of v implemented."
1529          def getV(self, p, v0):
1530           """           """
1531           Returns inner product of two elements v0 and v1 (overwrite).           return the value for v for a given p (overwrite)
1532
1533           @type v0: equal to the type of v           @param p: a pressure
1534           @type v1: equal to the type of v           @param v0: a initial guess for the value v to return.
1535           @rtype: equal to the type of v           @return: v given as M{v= A^{-1} (f-B^*p)}
1536           """           """
1537           raise NotImplementedError,"no inner product for v implemented."           raise NotImplementedError,"no v calculation implemented."
1538           pass
1539
1540          def Bv(self,v):
1541            """
1542            Returns Bv (overwrite).
1543
1544            @rtype: equal to the type of p
1545            @note: boundary conditions on p should be zero!
1546            """
1547            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_A(self,u,p):        def solve_AinvBt(self,p):
1559           """           """
1560           Solves M{Av=f-Au-B^*p} with accuracy L{self.getSubProblemTolerance()}           Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}
1561           (overwrite).           (overwrite).
1562
1563           @rtype: equal to the type of v           @param p: a pressure increment
1564             @return: the solution of M{Av=B^*p}
1565           @note: boundary conditions on v should be zero!           @note: boundary conditions on v should be zero!
1566           """           """
1567           raise NotImplementedError,"no operator A implemented."           raise NotImplementedError,"no operator A implemented."
1568
1569        def solve_prec(self,p):        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 (v,Bv) with v=A^-1B*p            dv=self.solve_AinvBt(p)
1587            #solve Av =B^*p as Av =f-Az-B^*(-p)            return ArithmeticTuple(dv,self.Bv(dv))
v=self.solve_A(self.__z,-p)
return ArithmeticTuple(v, self.B(v))
1588
1589        def __inner_PCG(self,p,r):        def __inner_PCG(self,p,r):
1590           return self.inner_pBv(p,r[1])           return self.inner_pBv(p,r[1])
# Line 1537  class HomogeneousSaddlePointProblem(obje Line 1592  class HomogeneousSaddlePointProblem(obje
1592        def __Msolve_PCG(self,r):        def __Msolve_PCG(self,r):
1593            return self.solve_prec(r[1])            return self.solve_prec(r[1])
1594        #=============================================================        #=============================================================
1595        def __Aprod_GMRES(self,x):        def __Aprod_GMRES(self,p):
1596            # return w,q  from v, p            return self.solve_prec(self.Bv(self.solve_AinvBt(p)))
1597            # solve Aw =Av+B^*p as Aw =f-A(z-v)-B^*(-p)        def __inner_GMRES(self,p0,p1):
1598            #  and  Sq=B(v-w)           return self.inner_p(p0,p1)
v=x[0]
p=x[1]
w=self.solve_A(self.__z-v,-p)
Bw=self.B(w-v)
q=self.solve_prec(Bw)
return ArithmeticTuple(w,q)

def __inner_GMRES(self,a1,a2):
return self.inner_p(a1[1],a2[1])+self.inner_v(a1[0],a2[0])
1599
1600        #=============================================================        #=============================================================
1601        def norm(self,v,p):        def norm_p(self,p):
1602          f=self.inner_p(p,p)+self.inner_v(v,v)            """
1603          if f<0:            calculates the norm of C{p}
1604              raise ValueError,"negative norm."
1605          return math.sqrt(f)            @param p: a pressure
1606              @return: the norm of C{p} using the inner product for pressure
1607        def solve(self,v,p,max_iter=20, verbose=False, show_details=False, useUzawa=True, iter_restart=20,max_correction_steps=3):            @rtype: C{float}
1608              """
1609              f=self.inner_p(p,p)
1610              if f<0: raise ValueError,"negative pressure norm."
1611              return math.sqrt(f)
1612          def adaptSubTolerance(self):
1613          """
1614          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 1566  class HomogeneousSaddlePointProblem(obje Line 1623  class HomogeneousSaddlePointProblem(obje
1623           @param p: initial guess for pressure           @param p: initial guess for pressure
1624           @type v: L{Data}           @type v: L{Data}
1625           @type p: L{Data}           @type p: L{Data}
1626           @param useUzawa: indicates the usage of the Uzawa scheme. Otherwise           @param usePCG: indicates the usage of the PCG rather than GMRES scheme.
the problem is solved in coupled form.
1627           @param max_iter: maximum number of iteration steps per correction           @param max_iter: maximum number of iteration steps per correction
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           @param max_correction_steps: maximum number of iteration steps in the           @type usePCG: C{bool}
attempt to get |Bv| to zero.
@return: new approximations for velocity and pressure
@type useUzawa: 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}
@type max_correction_steps: 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
#=====================================================================
# Az=f is solved as A(z-v)=f-Av-B^*0 (typically with z-v=0 on boundary)
self.__z=v+self.solve_A(v,p*0)
# tolerances:
1640           rtol=self.getTolerance()           rtol=self.getTolerance()
1641           atol=self.getAbsoluteTolerance()           atol=self.getAbsoluteTolerance()
1642           if useUzawa:       if self.adaptSubTolerance(): self.setSubProblemTolerance()
p0=self.solve_prec(self.B(self.__z))
f0=self.norm(self.__z,p0)
else:
f0=util.sqrt(self.inner_v(self.__z,self.__z))
if not f0>0: return self.__z, p*0
ATOL=rtol*f0+atol
if not ATOL>0: raise ValueError,"overall absolute tolerance needs to be positive."
if self.verbose: print "saddle point solver: initial residual %e, tolerance = %e."%(f0,ATOL)
# initialization
self.iter=0
1643           correction_step=0           correction_step=0
1644           converged=False           converged=False
1645           # initial guess:           while not converged:
1646           q=p*1                # calculate velocity for current pressure:
1647           u=v*1                v=self.getV(p,v)
1648           while not converged :                Bv=self.Bv(v)
1649              if useUzawa:                norm_v=self.norm_v(v)
1650                 # assume p is known: then v=z-A^-1B^*p                norm_Bv=self.norm_Bv(Bv)
1651                 #                ATOL=norm_v*rtol+atol
1652                 # which leads to BA^-1B^*p = Bz                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."
1654                 # note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv                if norm_Bv <= ATOL:
1655                 # we use the (v,Bv) to represent the residual                   converged=True
1656                 #                else:
1657                 # the norm of the right hand side Bv = f0                   correction_step+=1
1658                 #                   if correction_step>max_correction_steps:
1659                 #                  and the initial residual                        raise CorrectionFailed,"Given up after %d correction steps."%correction_step
1660                 #                   dp=self.solve_prec(Bv)
1661                 #     r=Bz-BA^-1B^*q = B(z-A^{-1}B^*q)=Bw with A(w-z)=Az-Az-B^*q = f -A*0 - B^{*}q                   if usePCG:
1662                 #                     norm2=self.inner_pBv(dp,Bv)
1663                 w=self.solve_A(self.__z,q)+self.__z                     if norm2<0: raise ValueError,"negative PCG norm."
1664                 if self.verbose: print "enter PCG method (iter_max=%s, atol=%e, subtolerance=%e)"%(max_iter,ATOL, self.getSubProblemTolerance())                     norm2=math.sqrt(norm2)
1665                 q,r=PCG(ArithmeticTuple(w,self.B(w)),self.__Aprod_PCG,q,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)                   else:
1666                 u=r[0]                     norm2=self.norm_p(dp)
1667                 Bu=r[1]                   ATOL_ITER=ATOL/norm_Bv*norm2*0.5
1668              else:                   if self.verbose: print "HomogeneousSaddlePointProblem: tolerance for solver: %e"%ATOL_ITER
1669                 #                   if usePCG:
1670                 #  with v=dv+z                         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:
1672                 #   A v + B p = f                         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                 #   B v       = 0           if self.verbose: print "HomogeneousSaddlePointProblem: tolerance reached."
1674                 #       return v,p
# apply the preconditioner [[A^{-1} 0][(S^{-1} B A^{-1})  -S^{-1}]]
#
w=self.solve_A(u,q)
if self.verbose: print "enter GMRES (iter_max=%s, atol=%e, subtolerance=%e)"%(max_iter,ATOL,self.getSubProblemTolerance())
x=GMRES(ArithmeticTuple(w,self.solve_prec(self.B(u+w))),self.__Aprod_GMRES, ArithmeticTuple(u,q), \
self.__inner_GMRES,atol=ATOL, rtol=0.,iter_max=max_iter, iter_restart=iter_restart, verbose=self.verbose)
u=x[0]
q=x[1]
Bu=self.B(u)
# now we check if |Bu| ~ 0 or more precise |Bu|_p  <= rtol * |v|_v
nu=self.inner_v(u,u)
p2=self.solve_prec(Bu)
nBu=self.inner_p(p2,p2)
if not nu>=0 and not nBu>=0: raise NegativeNorm,"negative norm."
nu=math.sqrt(nu)
nBu=math.sqrt(nBu)
if self.verbose: print "saddle point solver: norm v= %e (Bv = %e)"%(nu,nBu)
QTOL=atol+nu*rtol
if nBu <= QTOL:
converged=True
else:
ATOL=QTOL/nBu*ATOL*0.3
if self.verbose: print "correction step %s: tolerance reduced to %e."%(correction_step,ATOL)
converged=False
correction_step+=1
if correction_step>max_correction_steps:
raise CorrectionFailed,"Given up after %d correction steps."%correction_step
if self.verbose: print "saddle point solver: tolerance reached."
return u,q
1675
1676        #========================================================================        #========================================================================
1677        def setTolerance(self,tolerance=1.e-4):        def setTolerance(self,tolerance=1.e-4):
# Line 1678  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 1709  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 1752  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.2245 changed lines Added in v.2548

 ViewVC Help Powered by ViewVC 1.1.26