/[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 2793 by gross, Tue Dec 1 06:10:10 2009 UTC revision 3675 by jfenwick, Thu Nov 17 00:53:38 2011 UTC
# Line 1  Line 1 
1    
2  ########################################################  ########################################################
3  #  #
4  # Copyright (c) 2003-2009 by University of Queensland  # Copyright (c) 2003-2010 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 11  Line 11 
11  #  #
12  ########################################################  ########################################################
13    
14  __copyright__="""Copyright (c) 2003-2009 by University of Queensland  __copyright__="""Copyright (c) 2003-2010 by University of Queensland
15  Earth Systems Science Computational Center (ESSCC)  Earth Systems Science Computational Center (ESSCC)
16  http://www.uq.edu.au/esscc  http://www.uq.edu.au/esscc
17  Primary Business: Queensland, Australia"""  Primary Business: Queensland, Australia"""
# Line 153  class Projector: Line 153  class Projector:
153      """      """
154      return self.__pde.getSolverOptions()      return self.__pde.getSolverOptions()
155    
156      def getValue(self, input_data):
157        """
158        Projects ``input_data`` onto a continuous function.
159    
160        :param input_data: the data to be projected
161        """
162        return self(input_data)
163    
164    def __call__(self, input_data):    def __call__(self, input_data):
165      """      """
166      Projects ``input_data`` onto a continuous function.      Projects ``input_data`` onto a continuous function.
# Line 332  class Locator: Line 340  class Locator:
340         iterative=False         iterative=False
341         if isinstance(x, list):         if isinstance(x, list):
342             if len(x)==0:             if len(x)==0:
343                raise "ValueError", "At least one point must be given."                raise ValueError, "At least one point must be given."
344             try:             try:
345               iter(x[0])               iter(x[0])
346               iterative=True               iterative=True
347             except TypeError:             except TypeError:
348               iterative=False               iterative=False
349           xxx=self.__function_space.getX()
350         if iterative:         if iterative:
351             self.__id=[]             self.__id=[]
352             for p in x:             for p in x:
353                self.__id.append(util.length(self.__function_space.getX()-p[:self.__function_space.getDim()]).minGlobalDataPoint())                self.__id.append(util.length(xxx-p[:self.__function_space.getDim()]).minGlobalDataPoint())
354         else:         else:
355             self.__id=util.length(self.__function_space.getX()-x[:self.__function_space.getDim()]).minGlobalDataPoint()             self.__id=util.length(xxx-x[:self.__function_space.getDim()]).minGlobalDataPoint()
356    
357       def __str__(self):       def __str__(self):
358         """         """
# Line 400  class Locator: Line 409  class Locator:
409          Returns the value of ``data`` at the Locator if ``data`` is a `Data`          Returns the value of ``data`` at the Locator if ``data`` is a `Data`
410          object otherwise the object is returned.          object otherwise the object is returned.
411          """          """
412            print "Norris"
413          if isinstance(data,escript.Data):          if isinstance(data,escript.Data):
414             dat=util.interpolate(data,self.getFunctionSpace())             dat=util.interpolate(data,self.getFunctionSpace())
415             id=self.getId()             id=self.getId()
416             r=data.getRank()             r=data.getRank()
417             if isinstance(id,list):             if isinstance(id,list):
418                 out=[]                 out=[]
419                   print "Zorro"
420                 for i in id:                 for i in id:
421                    o=numpy.array(dat.getTupleForGlobalDataPoint(*i))                    o=numpy.array(dat.getTupleForGlobalDataPoint(*i))
422                    if data.getRank()==0:                    if data.getRank()==0:
# Line 414  class Locator: Line 425  class Locator:
425                       out.append(o)                       out.append(o)
426                 return out                 return out
427             else:             else:
428             print "Thom"
429               out=numpy.array(dat.getTupleForGlobalDataPoint(*id))               out=numpy.array(dat.getTupleForGlobalDataPoint(*id))
430               if data.getRank()==0:               if data.getRank()==0:
431                  return out[0]                  return out[0]
# Line 421  class Locator: Line 433  class Locator:
433                  return out                  return out
434          else:          else:
435             return data             return data
436              
437    #     def setValue(self, data, v):
438    #       """
439    #       Sets the value of the ``data`` at the Locator.
440    #       """
441    #       data.expand()   # Need to ensure that this is done globally
442    #       if isinstance(data, escript.Data):
443    #    id=self.getId()
444    #    if isinstance(id, list):
445    #      for i in id:
446    #         data._setTupleForGlobalDataPoint(i[1], i[0], v)
447    #    else:
448    #      data._setTupleForGlobalDataPoint(id[1], id[0], v)
449    #       else:
450    #    raise TypeError, "setValue: Invalid argument type."
451    
452    
453  def getInfLocator(arg):  def getInfLocator(arg):
# Line 1507  class HomogeneousSaddlePointProblem(obje Line 1534  class HomogeneousSaddlePointProblem(obje
1534          self.resetControlParameters()          self.resetControlParameters()
1535          self.setTolerance()          self.setTolerance()
1536          self.setAbsoluteTolerance()          self.setAbsoluteTolerance()
1537        def resetControlParameters(self,gamma=0.85, gamma_min=1.e-2,chi_max=0.1, omega_div=0.2, omega_conv=1.1, rtol_min=1.e-7, rtol_max=0.9, chi=1., C_p=1., C_v=1., safety_factor=0.3):        def resetControlParameters(self, K_p=1., K_v=1., rtol_max=0.01, rtol_min = 1.e-7, chi_max=0.5, reduction_factor=0.3, theta = 0.1):
1538           """           """
1539           sets a control parameter           sets a control parameter
1540    
1541           :param gamma: ``1/(1-gamma)`` controls the perturbation of the converegence rate due to termination errors in the subproblems.           :param K_p: initial value for constant to adjust pressure tolerance
1542           :type gamma: ``float``           :type K_p: ``float``
1543           :param gamma_min: minimum value for ``gamma``.           :param K_v: initial value for constant to adjust velocity tolerance
1544           :type gamma_min: ``float``           :type K_v: ``float``
          :param chi_max: maximum tolerable converegence rate.  
          :type chi_max: ``float``  
          :param omega_div: reduction fact for ``gamma`` if no convergence is detected.  
          :type omega_div: ``float``  
          :param omega_conv: raise fact for ``gamma`` if convergence is detected.  
          :type omega_conv: ``float``  
          :param rtol_min: minimum relative tolerance used to calculate presssure and velocity increment.  
          :type rtol_min: ``float``  
1545           :param rtol_max: maximuim relative tolerance used to calculate presssure and velocity increment.           :param rtol_max: maximuim relative tolerance used to calculate presssure and velocity increment.
1546           :type rtol_max: ``float``           :type rtol_max: ``float``
1547           :param chi: initial convergence measure.           :param chi_max: maximum tolerable converegence rate.
1548           :type chi: ``float``           :type chi_max: ``float``
1549           :param C_p: initial value for constant to adjust pressure tolerance           :param reduction_factor: reduction factor for adjustment factors.
1550           :type C_p: ``float``           :type reduction_factor: ``float``
          :param C_v: initial value for constant to adjust velocity tolerance  
          :type C_v: ``float``  
          :param safety_factor: safety factor for addjustment of pressure and velocity tolerance from stopping criteria  
          :type safety_factor: ``float``  
1551           """           """
1552           self.setControlParameter(gamma, gamma_min ,chi_max , omega_div , omega_conv, rtol_min , rtol_max, chi,C_p, C_v,safety_factor)           self.setControlParameter(K_p, K_v, rtol_max, rtol_min, chi_max, reduction_factor, theta)
1553    
1554        def setControlParameter(self,gamma=None, gamma_min=None ,chi_max=None, omega_div=None, omega_conv=None, rtol_min=None, rtol_max=None, chi=None, C_p=None, C_v=None, safety_factor=None):        def setControlParameter(self,K_p=None, K_v=None, rtol_max=None, rtol_min=None, chi_max=None, reduction_factor=None, theta=None):
1555           """           """
1556           sets a control parameter           sets a control parameter
1557    
1558           :param gamma: ``1/(1-gamma)`` controls the perturbation of the converegence rate due to termination errors in the subproblems.  
1559           :type gamma: ``float``           :param K_p: initial value for constant to adjust pressure tolerance
1560           :param gamma_min: minimum value for ``gamma``.           :type K_p: ``float``
1561           :type gamma_min: ``float``           :param K_v: initial value for constant to adjust velocity tolerance
1562           :param chi_max: maximum tolerable converegence rate.           :type K_v: ``float``
          :type chi_max: ``float``  
          :param omega_div: reduction fact for ``gamma`` if no convergence is detected.  
          :type omega_div: ``float``  
          :param omega_conv: raise fact for ``gamma`` if convergence is detected.  
          :type omega_conv: ``float``  
          :param rtol_min: minimum relative tolerance used to calculate presssure and velocity increment.  
          :type rtol_min: ``float``  
1563           :param rtol_max: maximuim relative tolerance used to calculate presssure and velocity increment.           :param rtol_max: maximuim relative tolerance used to calculate presssure and velocity increment.
1564           :type rtol_max: ``float``           :type rtol_max: ``float``
1565           :param chi: initial convergence measure.           :param chi_max: maximum tolerable converegence rate.
1566           :type chi: ``float``           :type chi_max: ``float``
1567           :param C_p: initial value for constant to adjust pressure tolerance           :type reduction_factor: ``float``
1568           :type C_p: ``float``           """
1569           :param C_v: initial value for constant to adjust velocity tolerance           if not K_p == None:
1570           :type C_v: ``float``              if K_p<1:
1571           :param safety_factor: safety factor for addjustment of pressure and velocity tolerance from stopping criteria                 raise ValueError,"K_p need to be greater or equal to 1."
          :type safety_factor: ``float``  
          """  
          if not gamma == None:  
             if gamma<=0 or gamma>=1:  
                raise ValueError,"Initial gamma needs to be positive and less than 1."  
          else:  
             gamma=self.__gamma  
   
          if not gamma_min == None:  
             if gamma_min<=0 or gamma_min>=1:  
                raise ValueError,"gamma_min needs to be positive and less than 1."  
          else:  
             gamma_min = self.__gamma_min  
   
          if not chi_max == None:  
             if chi_max<=0 or chi_max>=1:  
                raise ValueError,"chi_max needs to be positive and less than 1."  
          else:  
             chi_max = self.__chi_max  
   
          if not omega_div == None:  
             if omega_div<=0 or omega_div >=1:  
                raise ValueError,"omega_div needs to be positive and less than 1."  
          else:  
             omega_div=self.__omega_div  
   
          if not omega_conv == None:  
             if omega_conv<1:  
                raise ValueError,"omega_conv needs to be greater or equal to 1."  
1572           else:           else:
1573              omega_conv=self.__omega_conv              K_p=self.__K_p
1574    
1575           if not rtol_min == None:           if not K_v == None:
1576              if rtol_min<=0 or rtol_min>=1:              if K_v<1:
1577                 raise ValueError,"rtol_min needs to be positive and less than 1."                 raise ValueError,"K_v need to be greater or equal to 1."
1578           else:           else:
1579              rtol_min=self.__rtol_min              K_v=self.__K_v
1580    
1581           if not rtol_max == None:           if not rtol_max == None:
1582              if rtol_max<=0 or rtol_max>=1:              if rtol_max<=0 or rtol_max>=1:
# Line 1605  class HomogeneousSaddlePointProblem(obje Line 1584  class HomogeneousSaddlePointProblem(obje
1584           else:           else:
1585              rtol_max=self.__rtol_max              rtol_max=self.__rtol_max
1586    
1587           if not chi == None:           if not rtol_min == None:
1588              if chi<=0 or chi>1:              if rtol_min<=0 or rtol_min>=1:
1589                 raise ValueError,"chi needs to be positive and less than 1."                 raise ValueError,"rtol_min needs to be positive and less than 1."
1590           else:           else:
1591              chi=self.__chi              rtol_min=self.__rtol_min
1592    
1593           if not C_p == None:           if not chi_max == None:
1594              if C_p<1:              if chi_max<=0 or chi_max>=1:
1595                 raise ValueError,"C_p need to be greater or equal to 1."                 raise ValueError,"chi_max needs to be positive and less than 1."
1596           else:           else:
1597              C_p=self.__C_p              chi_max = self.__chi_max
1598    
1599           if not C_v == None:           if not reduction_factor == None:
1600              if C_v<1:              if reduction_factor<=0 or reduction_factor>1:
1601                 raise ValueError,"C_v need to be greater or equal to 1."                 raise ValueError,"reduction_factor need to be between zero and one."
1602           else:           else:
1603              C_v=self.__C_v              reduction_factor=self.__reduction_factor
1604    
1605           if not safety_factor == None:           if not theta == None:
1606              if safety_factor<=0 or safety_factor>1:              if theta<=0 or theta>1:
1607                 raise ValueError,"safety_factor need to be between zero and one."                 raise ValueError,"theta need to be between zero and one."
1608           else:           else:
1609              safety_factor=self.__safety_factor              theta=self.__theta
1610    
1611           if gamma<gamma_min:           if rtol_min>=rtol_max:
1612                 raise ValueError,"gamma = %e needs to be greater or equal gamma_min = %e."%(gamma,gamma_min)               raise ValueError,"rtol_max = %e needs to be greater than rtol_min = %e"%(rtol_max,rtol_min)
          if rtol_max<=rtol_min:  
                raise ValueError,"rtol_max = %e needs to be greater rtol_min = %e."%(rtol_max,rtol_min)  
               
          self.__gamma = gamma  
          self.__gamma_min = gamma_min  
1613           self.__chi_max = chi_max           self.__chi_max = chi_max
          self.__omega_div = omega_div  
          self.__omega_conv = omega_conv  
          self.__rtol_min = rtol_min  
1614           self.__rtol_max = rtol_max           self.__rtol_max = rtol_max
1615           self.__chi = chi           self.__K_p = K_p
1616           self.__C_p = C_p           self.__K_v = K_v
1617           self.__C_v = C_v           self.__reduction_factor = reduction_factor
1618           self.__safety_factor = safety_factor           self.__theta = theta
1619             self.__rtol_min=rtol_min
1620    
1621        #=============================================================        #=============================================================
1622        def inner_pBv(self,p,Bv):        def inner_pBv(self,p,Bv):
# Line 1787  class HomogeneousSaddlePointProblem(obje Line 1759  class HomogeneousSaddlePointProblem(obje
1759           self.verbose=verbose           self.verbose=verbose
1760           rtol=self.getTolerance()           rtol=self.getTolerance()
1761           atol=self.getAbsoluteTolerance()           atol=self.getAbsoluteTolerance()
1762    
1763             K_p=self.__K_p
1764             K_v=self.__K_v
1765           correction_step=0           correction_step=0
1766           converged=False           converged=False
          error=None  
1767           chi=None           chi=None
1768           gamma=self.__gamma           eps=None
1769           C_p=self.__C_p  
1770           C_v=self.__C_v           if self.verbose: print "HomogeneousSaddlePointProblem: start iteration: rtol= %e, atol=%e"%(rtol, atol)
1771           while not converged:           while not converged:
1772                if error== None or chi == None:  
1773                    gamma_new=gamma/self.__omega_conv               # get tolerance for velecity increment:
1774                else:               if chi == None:
1775                   if chi < self.__chi_max:                  rtol_v=self.__rtol_max
1776                      gamma_new=min(max(gamma*self.__omega_conv,1-chi*error/(self.__safety_factor*ATOL)), 1-chi/self.__chi_max)               else:
1777                   else:                  rtol_v=min(chi/K_v,self.__rtol_max)
1778                      gamma_new=gamma*self.__omega_div               rtol_v=max(rtol_v, self.__rtol_min)
1779                gamma=max(gamma_new, self.__gamma_min)               if self.verbose: print "HomogeneousSaddlePointProblem: step %s: rtol_v= %e"%(correction_step,rtol_v)
1780                # calculate velocity for current pressure:               # get velocity increment:
1781                rtol_v=min(max(gamma/(1.+gamma)/C_v,self.__rtol_min), self.__rtol_max)               dv1=self.getDV(p,v,rtol_v)
1782                rtol_p=min(max(gamma/C_p, self.__rtol_min), self.__rtol_max)               v1=v+dv1
1783                self.__subtol=rtol_p**2               Bv1=self.Bv(v1, rtol_v)
1784                if self.verbose: print "HomogeneousSaddlePointProblem: step %s: gamma = %e, rtol_v= %e, rtol_p=%e"%(correction_step,gamma,rtol_v,rtol_p)               norm_Bv1=self.norm_Bv(Bv1)
1785                if self.verbose: print "HomogeneousSaddlePointProblem: subtolerance: %e"%self.__subtol               norm_dv1=self.norm_v(dv1)
1786                # calculate velocity for current pressure: A*dv= F-A*v-B*p               if self.verbose: print "HomogeneousSaddlePointProblem: step %s: norm_Bv1 = %e, norm_dv1 = %e"%(correction_step, norm_Bv1, norm_dv1)
1787                dv1=self.getDV(p,v,rtol_v)               if norm_dv1*self.__theta < norm_Bv1:
1788                v1=v+dv1                  # get tolerance for pressure increment:
1789                Bv1=self.Bv(v1, self.__subtol)                  large_Bv1=True
1790                norm_Bv1=self.norm_Bv(Bv1)                  if chi == None or eps == None:
1791                norm_dv1=self.norm_v(dv1)                     rtol_p=self.__rtol_max
1792                norm_v1=self.norm_v(v1)                  else:
1793                ATOL=norm_v1*rtol+atol                     rtol_p=min(chi**2*eps/K_p/norm_Bv1, self.__rtol_max)
1794                if self.verbose: print "HomogeneousSaddlePointProblem: step %s: Bv = %e, dv = %e, v=%e"%(correction_step,norm_Bv1, norm_dv1, norm_v1)                  self.__subtol=max(rtol_p**2, self.__rtol_min)
1795                if not ATOL>0: raise ValueError,"overall absolute tolerance needs to be positive."                  if self.verbose: print "HomogeneousSaddlePointProblem: step %s: rtol_p= %e"%(correction_step,rtol_p)
1796                if max(norm_Bv1, norm_dv1) <= ATOL:                  # now we solve for the pressure increment dp from B*A^{-1}B^* dp = Bv1
1797                    converged=True                  if usePCG:
                   v=v1  
               else:  
                   # now we solve for the pressure increment dp from B*A^{-1}B^* dp = Bv1  
                   if usePCG:  
1798                      dp,r,a_norm=PCG(ArithmeticTuple(v1,Bv1),self.__Aprod_PCG,0*p,self.__Msolve_PCG,self.__inner_PCG,atol=0, rtol=rtol_p,iter_max=max_iter, verbose=self.verbose)                      dp,r,a_norm=PCG(ArithmeticTuple(v1,Bv1),self.__Aprod_PCG,0*p,self.__Msolve_PCG,self.__inner_PCG,atol=0, rtol=rtol_p,iter_max=max_iter, verbose=self.verbose)
1799                      v2=r[0]                      v2=r[0]
1800                      Bv2=r[1]                      Bv2=r[1]
1801                    else:                  else:
1802                        # don't use!!!!
1803                      dp=GMRES(self.solve_prec(Bv1,self.__subtol),self.__Aprod_GMRES, 0*p, self.__inner_GMRES,atol=0, rtol=rtol_p,iter_max=max_iter, iter_restart=iter_restart, verbose=self.verbose)                      dp=GMRES(self.solve_prec(Bv1,self.__subtol),self.__Aprod_GMRES, 0*p, self.__inner_GMRES,atol=0, rtol=rtol_p,iter_max=max_iter, iter_restart=iter_restart, verbose=self.verbose)
1804                      dv2=self.solve_AinvBt(dp, self.__subtol)                      dv2=self.solve_AinvBt(dp, self.__subtol)
1805                      v2=v1-dv2                      v2=v1-dv2
1806                      Bv2=self.Bv(v2, self.__subtol)                      Bv2=self.Bv(v2, self.__subtol)
1807                    #                  p2=p+dp
1808                    # convergence indicators:               else:
1809                    #                  large_Bv1=False
1810                    norm_v2=self.norm_v(v2)                  v2=v1
1811                    norm_dv2=self.norm_v(v2-v)                  p2=p
1812                    error_new=max(norm_dv2, norm_Bv1)               # update business:
1813                    norm_Bv2=self.norm_Bv(Bv2)               norm_dv2=self.norm_v(v2-v)
1814                    if self.verbose: print "HomogeneousSaddlePointProblem: step %s (part 2): Bv = %e, dv = %e, v=%e"%(correction_step,norm_Bv2, norm_dv2, norm_v2)               norm_v2=self.norm_v(v2)
1815                    if error !=None:               if self.verbose: print "HomogeneousSaddlePointProblem: step %s: v2 = %e, norm_dv2 = %e"%(correction_step, norm_v2, self.norm_v(v2-v))
1816                        chi_new=error_new/error               eps, eps_old = max(norm_Bv1, norm_dv2), eps
1817                        if self.verbose: print "HomogeneousSaddlePointProblem: step %s: convergence rate = %e, est. error = %e"%(correction_step,chi_new, error_new)               if eps_old == None:
1818                        if chi != None:                    chi, chi_old = None, chi
1819                            gamma0=max(gamma, 1-chi/chi_new)               else:
1820                            C_p*=gamma0/gamma                    chi, chi_old = min(eps/ eps_old, self.__chi_max), chi
1821                            C_v*=gamma0/gamma*(1+gamma)/(1+gamma0)               if eps != None:
1822                        chi=chi_new                   if chi !=None:
1823                    else:                      if self.verbose: print "HomogeneousSaddlePointProblem: step %s: convergence rate = %e, correction = %e"%(correction_step,chi, eps)
1824                        if self.verbose: print "HomogeneousSaddlePointProblem: step %s: est. error = %e"%(correction_step, error_new)                   else:
1825                        if self.verbose: print "HomogeneousSaddlePointProblem: step %s: correction = %e"%(correction_step, eps)
1826                    error = error_new               if eps <= rtol*norm_v2+atol :
1827                    correction_step+=1                   converged = True
1828                    if correction_step>max_correction_steps:               else:
1829                          raise CorrectionFailed,"Given up after %d correction steps."%correction_step                   if correction_step>=max_correction_steps:
1830                    v,p=v2,p+dp                        raise CorrectionFailed,"Given up after %d correction steps."%correction_step
1831                     if chi_old!=None:
1832                        K_p=max(1,self.__reduction_factor*K_p,(chi-chi_old)/chi_old**2*K_p)
1833                        K_v=max(1,self.__reduction_factor*K_v,(chi-chi_old)/chi_old**2*K_p)
1834                        if self.verbose: print "HomogeneousSaddlePointProblem: step %s: new adjustment factor K = %e"%(correction_step,K_p)
1835                 correction_step+=1
1836                 v,p =v2, p2
1837           if self.verbose: print "HomogeneousSaddlePointProblem: tolerance reached after %s steps."%correction_step           if self.verbose: print "HomogeneousSaddlePointProblem: tolerance reached after %s steps."%correction_step
1838       return v,p       return v,p
   
1839        #========================================================================        #========================================================================
1840        def setTolerance(self,tolerance=1.e-4):        def setTolerance(self,tolerance=1.e-4):
1841           """           """

Legend:
Removed from v.2793  
changed lines
  Added in v.3675

  ViewVC Help
Powered by ViewVC 1.1.26