/[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 3574 by jfenwick, Fri Sep 2 07:19:31 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 421  class Locator: Line 430  class Locator:
430                  return out                  return out
431          else:          else:
432             return data             return data
433              
434         def setValue(self, data, v):
435           """
436           Sets the value of the ``data`` at the Locator.
437           """
438           data.expand()    # Need to ensure that this is done globally
439           if isinstance(data, escript.Data):
440         id=self.getId()
441         if isinstance(id, list):
442           for i in id:
443              data._setTupleForGlobalDataPoint(i[1], i[0], v)
444         else:
445           data._setTupleForGlobalDataPoint(id[1], id[0], v)
446           else:
447         raise TypeError, "setValue: Invalid argument type."
448    
449    
450  def getInfLocator(arg):  def getInfLocator(arg):
# Line 1507  class HomogeneousSaddlePointProblem(obje Line 1531  class HomogeneousSaddlePointProblem(obje
1531          self.resetControlParameters()          self.resetControlParameters()
1532          self.setTolerance()          self.setTolerance()
1533          self.setAbsoluteTolerance()          self.setAbsoluteTolerance()
1534        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):
1535           """           """
1536           sets a control parameter           sets a control parameter
1537    
1538           :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
1539           :type gamma: ``float``           :type K_p: ``float``
1540           :param gamma_min: minimum value for ``gamma``.           :param K_v: initial value for constant to adjust velocity tolerance
1541           :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``  
1542           :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.
1543           :type rtol_max: ``float``           :type rtol_max: ``float``
1544           :param chi: initial convergence measure.           :param chi_max: maximum tolerable converegence rate.
1545           :type chi: ``float``           :type chi_max: ``float``
1546           :param C_p: initial value for constant to adjust pressure tolerance           :param reduction_factor: reduction factor for adjustment factors.
1547           :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``  
1548           """           """
1549           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)
1550    
1551        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):
1552           """           """
1553           sets a control parameter           sets a control parameter
1554    
1555           :param gamma: ``1/(1-gamma)`` controls the perturbation of the converegence rate due to termination errors in the subproblems.  
1556           :type gamma: ``float``           :param K_p: initial value for constant to adjust pressure tolerance
1557           :param gamma_min: minimum value for ``gamma``.           :type K_p: ``float``
1558           :type gamma_min: ``float``           :param K_v: initial value for constant to adjust velocity tolerance
1559           :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``  
1560           :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.
1561           :type rtol_max: ``float``           :type rtol_max: ``float``
1562           :param chi: initial convergence measure.           :param chi_max: maximum tolerable converegence rate.
1563           :type chi: ``float``           :type chi_max: ``float``
1564           :param C_p: initial value for constant to adjust pressure tolerance           :type reduction_factor: ``float``
1565           :type C_p: ``float``           """
1566           :param C_v: initial value for constant to adjust velocity tolerance           if not K_p == None:
1567           :type C_v: ``float``              if K_p<1:
1568           :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."  
1569           else:           else:
1570              omega_conv=self.__omega_conv              K_p=self.__K_p
1571    
1572           if not rtol_min == None:           if not K_v == None:
1573              if rtol_min<=0 or rtol_min>=1:              if K_v<1:
1574                 raise ValueError,"rtol_min needs to be positive and less than 1."                 raise ValueError,"K_v need to be greater or equal to 1."
1575           else:           else:
1576              rtol_min=self.__rtol_min              K_v=self.__K_v
1577    
1578           if not rtol_max == None:           if not rtol_max == None:
1579              if rtol_max<=0 or rtol_max>=1:              if rtol_max<=0 or rtol_max>=1:
# Line 1605  class HomogeneousSaddlePointProblem(obje Line 1581  class HomogeneousSaddlePointProblem(obje
1581           else:           else:
1582              rtol_max=self.__rtol_max              rtol_max=self.__rtol_max
1583    
1584           if not chi == None:           if not rtol_min == None:
1585              if chi<=0 or chi>1:              if rtol_min<=0 or rtol_min>=1:
1586                 raise ValueError,"chi needs to be positive and less than 1."                 raise ValueError,"rtol_min needs to be positive and less than 1."
1587           else:           else:
1588              chi=self.__chi              rtol_min=self.__rtol_min
1589    
1590           if not C_p == None:           if not chi_max == None:
1591              if C_p<1:              if chi_max<=0 or chi_max>=1:
1592                 raise ValueError,"C_p need to be greater or equal to 1."                 raise ValueError,"chi_max needs to be positive and less than 1."
1593           else:           else:
1594              C_p=self.__C_p              chi_max = self.__chi_max
1595    
1596           if not C_v == None:           if not reduction_factor == None:
1597              if C_v<1:              if reduction_factor<=0 or reduction_factor>1:
1598                 raise ValueError,"C_v need to be greater or equal to 1."                 raise ValueError,"reduction_factor need to be between zero and one."
1599           else:           else:
1600              C_v=self.__C_v              reduction_factor=self.__reduction_factor
1601    
1602           if not safety_factor == None:           if not theta == None:
1603              if safety_factor<=0 or safety_factor>1:              if theta<=0 or theta>1:
1604                 raise ValueError,"safety_factor need to be between zero and one."                 raise ValueError,"theta need to be between zero and one."
1605           else:           else:
1606              safety_factor=self.__safety_factor              theta=self.__theta
1607    
1608           if gamma<gamma_min:           if rtol_min>=rtol_max:
1609                 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  
1610           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  
1611           self.__rtol_max = rtol_max           self.__rtol_max = rtol_max
1612           self.__chi = chi           self.__K_p = K_p
1613           self.__C_p = C_p           self.__K_v = K_v
1614           self.__C_v = C_v           self.__reduction_factor = reduction_factor
1615           self.__safety_factor = safety_factor           self.__theta = theta
1616             self.__rtol_min=rtol_min
1617    
1618        #=============================================================        #=============================================================
1619        def inner_pBv(self,p,Bv):        def inner_pBv(self,p,Bv):
# Line 1787  class HomogeneousSaddlePointProblem(obje Line 1756  class HomogeneousSaddlePointProblem(obje
1756           self.verbose=verbose           self.verbose=verbose
1757           rtol=self.getTolerance()           rtol=self.getTolerance()
1758           atol=self.getAbsoluteTolerance()           atol=self.getAbsoluteTolerance()
1759    
1760             K_p=self.__K_p
1761             K_v=self.__K_v
1762           correction_step=0           correction_step=0
1763           converged=False           converged=False
          error=None  
1764           chi=None           chi=None
1765           gamma=self.__gamma           eps=None
1766           C_p=self.__C_p  
1767           C_v=self.__C_v           if self.verbose: print "HomogeneousSaddlePointProblem: start iteration: rtol= %e, atol=%e"%(rtol, atol)
1768           while not converged:           while not converged:
1769                if error== None or chi == None:  
1770                    gamma_new=gamma/self.__omega_conv               # get tolerance for velecity increment:
1771                else:               if chi == None:
1772                   if chi < self.__chi_max:                  rtol_v=self.__rtol_max
1773                      gamma_new=min(max(gamma*self.__omega_conv,1-chi*error/(self.__safety_factor*ATOL)), 1-chi/self.__chi_max)               else:
1774                   else:                  rtol_v=min(chi/K_v,self.__rtol_max)
1775                      gamma_new=gamma*self.__omega_div               rtol_v=max(rtol_v, self.__rtol_min)
1776                gamma=max(gamma_new, self.__gamma_min)               if self.verbose: print "HomogeneousSaddlePointProblem: step %s: rtol_v= %e"%(correction_step,rtol_v)
1777                # calculate velocity for current pressure:               # get velocity increment:
1778                rtol_v=min(max(gamma/(1.+gamma)/C_v,self.__rtol_min), self.__rtol_max)               dv1=self.getDV(p,v,rtol_v)
1779                rtol_p=min(max(gamma/C_p, self.__rtol_min), self.__rtol_max)               v1=v+dv1
1780                self.__subtol=rtol_p**2               Bv1=self.Bv(v1, rtol_v)
1781                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)
1782                if self.verbose: print "HomogeneousSaddlePointProblem: subtolerance: %e"%self.__subtol               norm_dv1=self.norm_v(dv1)
1783                # 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)
1784                dv1=self.getDV(p,v,rtol_v)               if norm_dv1*self.__theta < norm_Bv1:
1785                v1=v+dv1                  # get tolerance for pressure increment:
1786                Bv1=self.Bv(v1, self.__subtol)                  large_Bv1=True
1787                norm_Bv1=self.norm_Bv(Bv1)                  if chi == None or eps == None:
1788                norm_dv1=self.norm_v(dv1)                     rtol_p=self.__rtol_max
1789                norm_v1=self.norm_v(v1)                  else:
1790                ATOL=norm_v1*rtol+atol                     rtol_p=min(chi**2*eps/K_p/norm_Bv1, self.__rtol_max)
1791                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)
1792                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)
1793                if max(norm_Bv1, norm_dv1) <= ATOL:                  # now we solve for the pressure increment dp from B*A^{-1}B^* dp = Bv1
1794                    converged=True                  if usePCG:
                   v=v1  
               else:  
                   # now we solve for the pressure increment dp from B*A^{-1}B^* dp = Bv1  
                   if usePCG:  
1795                      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)
1796                      v2=r[0]                      v2=r[0]
1797                      Bv2=r[1]                      Bv2=r[1]
1798                    else:                  else:
1799                        # don't use!!!!
1800                      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)
1801                      dv2=self.solve_AinvBt(dp, self.__subtol)                      dv2=self.solve_AinvBt(dp, self.__subtol)
1802                      v2=v1-dv2                      v2=v1-dv2
1803                      Bv2=self.Bv(v2, self.__subtol)                      Bv2=self.Bv(v2, self.__subtol)
1804                    #                  p2=p+dp
1805                    # convergence indicators:               else:
1806                    #                  large_Bv1=False
1807                    norm_v2=self.norm_v(v2)                  v2=v1
1808                    norm_dv2=self.norm_v(v2-v)                  p2=p
1809                    error_new=max(norm_dv2, norm_Bv1)               # update business:
1810                    norm_Bv2=self.norm_Bv(Bv2)               norm_dv2=self.norm_v(v2-v)
1811                    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)
1812                    if error !=None:               if self.verbose: print "HomogeneousSaddlePointProblem: step %s: v2 = %e, norm_dv2 = %e"%(correction_step, norm_v2, self.norm_v(v2-v))
1813                        chi_new=error_new/error               eps, eps_old = max(norm_Bv1, norm_dv2), eps
1814                        if self.verbose: print "HomogeneousSaddlePointProblem: step %s: convergence rate = %e, est. error = %e"%(correction_step,chi_new, error_new)               if eps_old == None:
1815                        if chi != None:                    chi, chi_old = None, chi
1816                            gamma0=max(gamma, 1-chi/chi_new)               else:
1817                            C_p*=gamma0/gamma                    chi, chi_old = min(eps/ eps_old, self.__chi_max), chi
1818                            C_v*=gamma0/gamma*(1+gamma)/(1+gamma0)               if eps != None:
1819                        chi=chi_new                   if chi !=None:
1820                    else:                      if self.verbose: print "HomogeneousSaddlePointProblem: step %s: convergence rate = %e, correction = %e"%(correction_step,chi, eps)
1821                        if self.verbose: print "HomogeneousSaddlePointProblem: step %s: est. error = %e"%(correction_step, error_new)                   else:
1822                        if self.verbose: print "HomogeneousSaddlePointProblem: step %s: correction = %e"%(correction_step, eps)
1823                    error = error_new               if eps <= rtol*norm_v2+atol :
1824                    correction_step+=1                   converged = True
1825                    if correction_step>max_correction_steps:               else:
1826                          raise CorrectionFailed,"Given up after %d correction steps."%correction_step                   if correction_step>=max_correction_steps:
1827                    v,p=v2,p+dp                        raise CorrectionFailed,"Given up after %d correction steps."%correction_step
1828                     if chi_old!=None:
1829                        K_p=max(1,self.__reduction_factor*K_p,(chi-chi_old)/chi_old**2*K_p)
1830                        K_v=max(1,self.__reduction_factor*K_v,(chi-chi_old)/chi_old**2*K_p)
1831                        if self.verbose: print "HomogeneousSaddlePointProblem: step %s: new adjustment factor K = %e"%(correction_step,K_p)
1832                 correction_step+=1
1833                 v,p =v2, p2
1834           if self.verbose: print "HomogeneousSaddlePointProblem: tolerance reached after %s steps."%correction_step           if self.verbose: print "HomogeneousSaddlePointProblem: tolerance reached after %s steps."%correction_step
1835       return v,p       return v,p
   
1836        #========================================================================        #========================================================================
1837        def setTolerance(self,tolerance=1.e-4):        def setTolerance(self,tolerance=1.e-4):
1838           """           """

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

  ViewVC Help
Powered by ViewVC 1.1.26