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

Legend:
Removed from v.2842  
changed lines
  Added in v.2843

  ViewVC Help
Powered by ViewVC 1.1.26