/[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 2881 by jfenwick, Thu Jan 28 02:03:15 2010 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 1507  class HomogeneousSaddlePointProblem(obje Line 1515  class HomogeneousSaddlePointProblem(obje
1515          self.resetControlParameters()          self.resetControlParameters()
1516          self.setTolerance()          self.setTolerance()
1517          self.setAbsoluteTolerance()          self.setAbsoluteTolerance()
1518        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):
1519           """           """
1520           sets a control parameter           sets a control parameter
1521    
1522           :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
1523           :type gamma: ``float``           :type K_p: ``float``
1524           :param gamma_min: minimum value for ``gamma``.           :param K_v: initial value for constant to adjust velocity tolerance
1525           :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``  
1526           :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.
1527           :type rtol_max: ``float``           :type rtol_max: ``float``
1528           :param chi: initial convergence measure.           :param chi_max: maximum tolerable converegence rate.
1529           :type chi: ``float``           :type chi_max: ``float``
1530           :param C_p: initial value for constant to adjust pressure tolerance           :param reduction_factor: reduction factor for adjustment factors.
1531           :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``  
1532           """           """
1533           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)
1534    
1535        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):
1536           """           """
1537           sets a control parameter           sets a control parameter
1538    
1539           :param gamma: ``1/(1-gamma)`` controls the perturbation of the converegence rate due to termination errors in the subproblems.  
1540           :type gamma: ``float``           :param K_p: initial value for constant to adjust pressure tolerance
1541           :param gamma_min: minimum value for ``gamma``.           :type K_p: ``float``
1542           :type gamma_min: ``float``           :param K_v: initial value for constant to adjust velocity tolerance
1543           :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``  
1544           :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.
1545           :type rtol_max: ``float``           :type rtol_max: ``float``
1546           :param chi: initial convergence measure.           :param chi_max: maximum tolerable converegence rate.
1547           :type chi: ``float``           :type chi_max: ``float``
1548           :param C_p: initial value for constant to adjust pressure tolerance           :type reduction_factor: ``float``
1549           :type C_p: ``float``           """
1550           :param C_v: initial value for constant to adjust velocity tolerance           if not K_p == None:
1551           :type C_v: ``float``              if K_p<1:
1552           :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."  
1553           else:           else:
1554              omega_conv=self.__omega_conv              K_p=self.__K_p
1555    
1556           if not rtol_min == None:           if not K_v == None:
1557              if rtol_min<=0 or rtol_min>=1:              if K_v<1:
1558                 raise ValueError,"rtol_min needs to be positive and less than 1."                 raise ValueError,"K_v need to be greater or equal to 1."
1559           else:           else:
1560              rtol_min=self.__rtol_min              K_v=self.__K_v
1561    
1562           if not rtol_max == None:           if not rtol_max == None:
1563              if rtol_max<=0 or rtol_max>=1:              if rtol_max<=0 or rtol_max>=1:
# Line 1605  class HomogeneousSaddlePointProblem(obje Line 1565  class HomogeneousSaddlePointProblem(obje
1565           else:           else:
1566              rtol_max=self.__rtol_max              rtol_max=self.__rtol_max
1567    
1568           if not chi == None:           if not rtol_min == None:
1569              if chi<=0 or chi>1:              if rtol_min<=0 or rtol_min>=1:
1570                 raise ValueError,"chi needs to be positive and less than 1."                 raise ValueError,"rtol_min needs to be positive and less than 1."
1571           else:           else:
1572              chi=self.__chi              rtol_min=self.__rtol_min
1573    
1574           if not C_p == None:           if not chi_max == None:
1575              if C_p<1:              if chi_max<=0 or chi_max>=1:
1576                 raise ValueError,"C_p need to be greater or equal to 1."                 raise ValueError,"chi_max needs to be positive and less than 1."
1577           else:           else:
1578              C_p=self.__C_p              chi_max = self.__chi_max
1579    
1580           if not C_v == None:           if not reduction_factor == None:
1581              if C_v<1:              if reduction_factor<=0 or reduction_factor>1:
1582                 raise ValueError,"C_v need to be greater or equal to 1."                 raise ValueError,"reduction_factor need to be between zero and one."
1583           else:           else:
1584              C_v=self.__C_v              reduction_factor=self.__reduction_factor
1585    
1586           if not safety_factor == None:           if not theta == None:
1587              if safety_factor<=0 or safety_factor>1:              if theta<=0 or theta>1:
1588                 raise ValueError,"safety_factor need to be between zero and one."                 raise ValueError,"theta need to be between zero and one."
1589           else:           else:
1590              safety_factor=self.__safety_factor              theta=self.__theta
1591    
1592           if gamma<gamma_min:           if rtol_min>=rtol_max:
1593                 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  
1594           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  
1595           self.__rtol_max = rtol_max           self.__rtol_max = rtol_max
1596           self.__chi = chi           self.__K_p = K_p
1597           self.__C_p = C_p           self.__K_v = K_v
1598           self.__C_v = C_v           self.__reduction_factor = reduction_factor
1599           self.__safety_factor = safety_factor           self.__theta = theta
1600             self.__rtol_min=rtol_min
1601    
1602        #=============================================================        #=============================================================
1603        def inner_pBv(self,p,Bv):        def inner_pBv(self,p,Bv):
# Line 1787  class HomogeneousSaddlePointProblem(obje Line 1740  class HomogeneousSaddlePointProblem(obje
1740           self.verbose=verbose           self.verbose=verbose
1741           rtol=self.getTolerance()           rtol=self.getTolerance()
1742           atol=self.getAbsoluteTolerance()           atol=self.getAbsoluteTolerance()
1743    
1744             K_p=self.__K_p
1745             K_v=self.__K_v
1746           correction_step=0           correction_step=0
1747           converged=False           converged=False
          error=None  
1748           chi=None           chi=None
1749           gamma=self.__gamma           eps=None
1750           C_p=self.__C_p  
1751           C_v=self.__C_v           if self.verbose: print "HomogeneousSaddlePointProblem: start iteration: rtol= %e, atol=%e"%(rtol, atol)
1752           while not converged:           while not converged:
1753                if error== None or chi == None:  
1754                    gamma_new=gamma/self.__omega_conv               # get tolerance for velecity increment:
1755                else:               if chi == None:
1756                   if chi < self.__chi_max:                  rtol_v=self.__rtol_max
1757                      gamma_new=min(max(gamma*self.__omega_conv,1-chi*error/(self.__safety_factor*ATOL)), 1-chi/self.__chi_max)               else:
1758                   else:                  rtol_v=min(chi/K_v,self.__rtol_max)
1759                      gamma_new=gamma*self.__omega_div               rtol_v=max(rtol_v, self.__rtol_min)
1760                gamma=max(gamma_new, self.__gamma_min)               if self.verbose: print "HomogeneousSaddlePointProblem: step %s: rtol_v= %e"%(correction_step,rtol_v)
1761                # calculate velocity for current pressure:               # get velocity increment:
1762                rtol_v=min(max(gamma/(1.+gamma)/C_v,self.__rtol_min), self.__rtol_max)               dv1=self.getDV(p,v,rtol_v)
1763                rtol_p=min(max(gamma/C_p, self.__rtol_min), self.__rtol_max)               v1=v+dv1
1764                self.__subtol=rtol_p**2               Bv1=self.Bv(v1, rtol_v)
1765                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)
1766                if self.verbose: print "HomogeneousSaddlePointProblem: subtolerance: %e"%self.__subtol               norm_dv1=self.norm_v(dv1)
1767                # 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)
1768                dv1=self.getDV(p,v,rtol_v)               if norm_dv1*self.__theta < norm_Bv1:
1769                v1=v+dv1                  # get tolerance for pressure increment:
1770                Bv1=self.Bv(v1, self.__subtol)                  large_Bv1=True
1771                norm_Bv1=self.norm_Bv(Bv1)                  if chi == None or eps == None:
1772                norm_dv1=self.norm_v(dv1)                     rtol_p=self.__rtol_max
1773                norm_v1=self.norm_v(v1)                  else:
1774                ATOL=norm_v1*rtol+atol                     rtol_p=min(chi**2*eps/K_p/norm_Bv1, self.__rtol_max)
1775                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)
1776                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)
1777                if max(norm_Bv1, norm_dv1) <= ATOL:                  # now we solve for the pressure increment dp from B*A^{-1}B^* dp = Bv1
1778                    converged=True                  if usePCG:
                   v=v1  
               else:  
                   # now we solve for the pressure increment dp from B*A^{-1}B^* dp = Bv1  
                   if usePCG:  
1779                      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)
1780                      v2=r[0]                      v2=r[0]
1781                      Bv2=r[1]                      Bv2=r[1]
1782                    else:                  else:
1783                        # don't use!!!!
1784                      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)
1785                      dv2=self.solve_AinvBt(dp, self.__subtol)                      dv2=self.solve_AinvBt(dp, self.__subtol)
1786                      v2=v1-dv2                      v2=v1-dv2
1787                      Bv2=self.Bv(v2, self.__subtol)                      Bv2=self.Bv(v2, self.__subtol)
1788                    #                  p2=p+dp
1789                    # convergence indicators:               else:
1790                    #                  large_Bv1=False
1791                    norm_v2=self.norm_v(v2)                  v2=v1
1792                    norm_dv2=self.norm_v(v2-v)                  p2=p
1793                    error_new=max(norm_dv2, norm_Bv1)               # update business:
1794                    norm_Bv2=self.norm_Bv(Bv2)               norm_dv2=self.norm_v(v2-v)
1795                    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)
1796                    if error !=None:               if self.verbose: print "HomogeneousSaddlePointProblem: step %s: v2 = %e, norm_dv2 = %e"%(correction_step, norm_v2, self.norm_v(v2-v))
1797                        chi_new=error_new/error               eps, eps_old = max(norm_Bv1, norm_dv2), eps
1798                        if self.verbose: print "HomogeneousSaddlePointProblem: step %s: convergence rate = %e, est. error = %e"%(correction_step,chi_new, error_new)               if eps_old == None:
1799                        if chi != None:                    chi, chi_old = None, chi
1800                            gamma0=max(gamma, 1-chi/chi_new)               else:
1801                            C_p*=gamma0/gamma                    chi, chi_old = min(eps/ eps_old, self.__chi_max), chi
1802                            C_v*=gamma0/gamma*(1+gamma)/(1+gamma0)               if eps != None:
1803                        chi=chi_new                   if chi !=None:
1804                    else:                      if self.verbose: print "HomogeneousSaddlePointProblem: step %s: convergence rate = %e, correction = %e"%(correction_step,chi, eps)
1805                        if self.verbose: print "HomogeneousSaddlePointProblem: step %s: est. error = %e"%(correction_step, error_new)                   else:
1806                        if self.verbose: print "HomogeneousSaddlePointProblem: step %s: correction = %e"%(correction_step, eps)
1807                    error = error_new               if eps <= rtol*norm_v2+atol :
1808                    correction_step+=1                   converged = True
1809                    if correction_step>max_correction_steps:               else:
1810                          raise CorrectionFailed,"Given up after %d correction steps."%correction_step                   if correction_step>=max_correction_steps:
1811                    v,p=v2,p+dp                        raise CorrectionFailed,"Given up after %d correction steps."%correction_step
1812                     if chi_old!=None:
1813                        K_p=max(1,self.__reduction_factor*K_p,(chi-chi_old)/chi_old**2*K_p)
1814                        K_v=max(1,self.__reduction_factor*K_v,(chi-chi_old)/chi_old**2*K_p)
1815                        if self.verbose: print "HomogeneousSaddlePointProblem: step %s: new adjustment factor K = %e"%(correction_step,K_p)
1816                 correction_step+=1
1817                 v,p =v2, p2
1818           if self.verbose: print "HomogeneousSaddlePointProblem: tolerance reached after %s steps."%correction_step           if self.verbose: print "HomogeneousSaddlePointProblem: tolerance reached after %s steps."%correction_step
1819       return v,p       return v,p
   
1820        #========================================================================        #========================================================================
1821        def setTolerance(self,tolerance=1.e-4):        def setTolerance(self,tolerance=1.e-4):
1822           """           """

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

  ViewVC Help
Powered by ViewVC 1.1.26