/[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 3851 by jfenwick, Wed Jan 18 02:30:48 2012 UTC revision 3852 by jfenwick, Thu Mar 1 05:34:52 2012 UTC
# Line 452  def getInfLocator(arg): Line 452  def getInfLocator(arg):
452      Return a Locator for a point with the inf value over all arg.      Return a Locator for a point with the inf value over all arg.
453      """      """
454      if not isinstance(arg, escript.Data):      if not isinstance(arg, escript.Data):
455      raise TypeError("getInfLocator: Unknown argument type.")         raise TypeError("getInfLocator: Unknown argument type.")
456      a_inf=util.inf(arg)      a_inf=util.inf(arg)
457      loc=util.length(arg-a_inf).minGlobalDataPoint() # This gives us the location but not coords      loc=util.length(arg-a_inf).minGlobalDataPoint() # This gives us the location but not coords
458      x=arg.getFunctionSpace().getX()      x=arg.getFunctionSpace().getX()
# Line 464  def getSupLocator(arg): Line 464  def getSupLocator(arg):
464      Return a Locator for a point with the sup value over all arg.      Return a Locator for a point with the sup value over all arg.
465      """      """
466      if not isinstance(arg, escript.Data):      if not isinstance(arg, escript.Data):
467      raise TypeError("getInfLocator: Unknown argument type.")         raise TypeError("getInfLocator: Unknown argument type.")
468      a_inf=util.sup(arg)      a_inf=util.sup(arg)
469      loc=util.length(arg-a_inf).minGlobalDataPoint() # This gives us the location but not coords      loc=util.length(arg-a_inf).minGlobalDataPoint() # This gives us the location but not coords
470      x=arg.getFunctionSpace().getX()      x=arg.getFunctionSpace().getX()
# Line 581  def PCG(r, Aprod, x, Msolve, bilinearfor Line 581  def PCG(r, Aprod, x, Msolve, bilinearfor
581         alpha = rhat_dot_r / bilinearform(d, q)         alpha = rhat_dot_r / bilinearform(d, q)
582         x += alpha * d         x += alpha * d
583         if isinstance(q,ArithmeticTuple):         if isinstance(q,ArithmeticTuple):
584         r += q * (-alpha)      # Doing it the other way calls the float64.__mul__ not AT.__rmul__            r += q * (-alpha)      # Doing it the other way calls the float64.__mul__ not AT.__rmul__
585         else:         else:
586             r += (-alpha) * q             r += (-alpha) * q
587         rhat=Msolve(r)         rhat=Msolve(r)
# Line 744  def NewtonGMRES(defect, x, iter_max=100, Line 744  def NewtonGMRES(defect, x, iter_max=100,
744          #   adjust subtol_          #   adjust subtol_
745          #          #
746              if iter > 1:              if iter > 1:
747             rat=fnrm/fnrmo                 rat=fnrm/fnrmo
748                 subtol_old=subtol                 subtol_old=subtol
749             subtol=gamma*rat**2                 subtol=gamma*rat**2
750             if gamma*subtol_old**2 > .1: subtol=max(subtol,gamma*subtol_old**2)                 if gamma*subtol_old**2 > .1: subtol=max(subtol,gamma*subtol_old**2)
751             subtol=max(min(subtol,subtol_max), .5*stop_tol/fnrm)                 subtol=max(min(subtol,subtol_max), .5*stop_tol/fnrm)
752          #          #
753          # calculate newton increment xc          # calculate newton increment xc
754              #     if iter_max in __FDGMRES is reached MaxIterReached is thrown              #     if iter_max in __FDGMRES is reached MaxIterReached is thrown
# Line 766  def NewtonGMRES(defect, x, iter_max=100, Line 766  def NewtonGMRES(defect, x, iter_max=100,
766              else:              else:
767                 iter+=sub_iter                 iter+=sub_iter
768              # ====              # ====
769          x+=xc              x+=xc
770              F=defect(x)              F=defect(x)
771          iter+=1              iter+=1
772              fnrmo, fnrm=fnrm, defect.norm(F)              fnrmo, fnrm=fnrm, defect.norm(F)
773              if verbose: print(("             step %s: residual %e."%(iter,fnrm)))              if verbose: print(("             step %s: residual %e."%(iter,fnrm)))
774     if verbose: print(("NewtonGMRES: completed after %s steps."%iter))     if verbose: print(("NewtonGMRES: completed after %s steps."%iter))
# Line 787  def __givapp(c,s,vin): Line 787  def __givapp(c,s,vin):
787      else:      else:
788          for i in range(len(c)):          for i in range(len(c)):
789              w1=c[i]*vrot[i]-s[i]*vrot[i+1]              w1=c[i]*vrot[i]-s[i]*vrot[i+1]
790          w2=s[i]*vrot[i]+c[i]*vrot[i+1]              w2=s[i]*vrot[i]+c[i]*vrot[i+1]
791              vrot[i]=w1              vrot[i]=w1
792              vrot[i+1]=w2              vrot[i+1]=w2
793      return vrot      return vrot
# Line 869  def __FDGMRES(F0, defect, x0, atol, iter Line 869  def __FDGMRES(F0, defect, x0, atol, iter
869            i=i-1            i=i-1
870       xhat=v[iter-1]*y[iter-1]       xhat=v[iter-1]*y[iter-1]
871       for i in range(iter-1):       for i in range(iter-1):
872      xhat += v[i]*y[i]         xhat += v[i]*y[i]
873     else :     else :
874        xhat=v[0] * 0        xhat=v[0] * 0
875    
# Line 971  def _GMRESm(r, Aprod, x, bilinearform, a Line 971  def _GMRESm(r, Aprod, x, bilinearform, a
971     if verbose: print(("GMRES: initial residual %e (absolute tolerance = %e)"%(rho,atol)))     if verbose: print(("GMRES: initial residual %e (absolute tolerance = %e)"%(rho,atol)))
972     while not (rho<=atol or iter==iter_restart):     while not (rho<=atol or iter==iter_restart):
973    
974      if iter  >= iter_max: raise MaxIterReached("maximum number of %s steps reached."%iter_max)          if iter  >= iter_max: raise MaxIterReached("maximum number of %s steps reached."%iter_max)
975    
976          if P_R!=None:          if P_R!=None:
977              p=Aprod(P_R(v[iter]))              p=Aprod(P_R(v[iter]))
978          else:          else:
979          p=Aprod(v[iter])              p=Aprod(v[iter])
980      v.append(p)          v.append(p)
981    
982      v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))          v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
983    
984  # Modified Gram-Schmidt  # Modified Gram-Schmidt
985      for j in range(iter+1):          for j in range(iter+1):
986        h[j,iter]=bilinearform(v[j],v[iter+1])            h[j,iter]=bilinearform(v[j],v[iter+1])
987        v[iter+1]-=h[j,iter]*v[j]            v[iter+1]-=h[j,iter]*v[j]
988    
989      h[iter+1,iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))          h[iter+1,iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
990      v_norm2=h[iter+1,iter]          v_norm2=h[iter+1,iter]
991    
992  # Reorthogonalize if needed  # Reorthogonalize if needed
993      if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)          if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
994       for j in range(iter+1):           for j in range(iter+1):
995          hr=bilinearform(v[j],v[iter+1])              hr=bilinearform(v[j],v[iter+1])
996              h[j,iter]=h[j,iter]+hr              h[j,iter]=h[j,iter]+hr
997              v[iter+1] -= hr*v[j]              v[iter+1] -= hr*v[j]
998    
999       v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))           v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
1000       h[iter+1,iter]=v_norm2           h[iter+1,iter]=v_norm2
1001    
1002  #   watch out for happy breakdown  #   watch out for happy breakdown
1003          if not v_norm2 == 0:          if not v_norm2 == 0:
1004           v[iter+1]=v[iter+1]/h[iter+1,iter]           v[iter+1]=v[iter+1]/h[iter+1,iter]
1005    
1006  #   Form and store the information for the new Givens rotation  #   Form and store the information for the new Givens rotation
1007      if iter > 0: h[:iter+1,iter]=__givapp(c[:iter],s[:iter],h[:iter+1,iter])          if iter > 0: h[:iter+1,iter]=__givapp(c[:iter],s[:iter],h[:iter+1,iter])
1008      mu=math.sqrt(h[iter,iter]*h[iter,iter]+h[iter+1,iter]*h[iter+1,iter])          mu=math.sqrt(h[iter,iter]*h[iter,iter]+h[iter+1,iter]*h[iter+1,iter])
1009    
1010      if mu!=0 :          if mu!=0 :
1011          c[iter]=h[iter,iter]/mu                  c[iter]=h[iter,iter]/mu
1012          s[iter]=-h[iter+1,iter]/mu                  s[iter]=-h[iter+1,iter]/mu
1013          h[iter,iter]=c[iter]*h[iter,iter]-s[iter]*h[iter+1,iter]                  h[iter,iter]=c[iter]*h[iter,iter]-s[iter]*h[iter+1,iter]
1014          h[iter+1,iter]=0.0                  h[iter+1,iter]=0.0
1015                  gg=__givapp(c[iter],s[iter],[g[iter],g[iter+1]])                  gg=__givapp(c[iter],s[iter],[g[iter],g[iter+1]])
1016                  g[iter]=gg[0]                  g[iter]=gg[0]
1017                  g[iter+1]=gg[1]                  g[iter+1]=gg[1]
# Line 1019  def _GMRESm(r, Aprod, x, bilinearform, a Line 1019  def _GMRESm(r, Aprod, x, bilinearform, a
1019    
1020          rho=abs(g[iter+1])          rho=abs(g[iter+1])
1021          if verbose: print(("GMRES: iteration step %s: residual %e"%(iter,rho)))          if verbose: print(("GMRES: iteration step %s: residual %e"%(iter,rho)))
1022      iter+=1          iter+=1
1023    
1024  # At this point either iter > iter_max or rho < tol.  # At this point either iter > iter_max or rho < tol.
1025  # It's time to compute x and leave.  # It's time to compute x and leave.
# Line 1035  def _GMRESm(r, Aprod, x, bilinearform, a Line 1035  def _GMRESm(r, Aprod, x, bilinearform, a
1035            i=i-1            i=i-1
1036       xhat=v[iter-1]*y[iter-1]       xhat=v[iter-1]*y[iter-1]
1037       for i in range(iter-1):       for i in range(iter-1):
1038      xhat += v[i]*y[i]         xhat += v[i]*y[i]
1039     else:     else:
1040       xhat=v[0] * 0       xhat=v[0] * 0
1041     if P_R!=None:     if P_R!=None:
# Line 1143  def MINRES(r, Aprod, x, Msolve, bilinear Line 1143  def MINRES(r, Aprod, x, Msolve, bilinear
1143      # --------------------------------------------------------------------      # --------------------------------------------------------------------
1144      while not rnorm<=atol+rtol*Anorm*ynorm:    #  checks ||r|| < (||A|| ||x||) * TOL      while not rnorm<=atol+rtol*Anorm*ynorm:    #  checks ||r|| < (||A|| ||x||) * TOL
1145    
1146      if iter  >= iter_max: raise MaxIterReached("maximum number of %s steps reached."%iter_max)          if iter  >= iter_max: raise MaxIterReached("maximum number of %s steps reached."%iter_max)
1147          iter    = iter  +  1          iter    = iter  +  1
1148    
1149          #-----------------------------------------------------------------          #-----------------------------------------------------------------
# Line 1525  class HomogeneousSaddlePointProblem(obje Line 1525  class HomogeneousSaddlePointProblem(obje
1525        *A* may depend weakly on *v* and *p*.        *A* may depend weakly on *v* and *p*.
1526        """        """
1527        def __init__(self, **kwargs):        def __init__(self, **kwargs):
1528      """          """
1529      initializes the saddle point problem          initializes the saddle point problem
1530      """          """
1531          self.resetControlParameters()          self.resetControlParameters()
1532          self.setTolerance()          self.setTolerance()
1533          self.setAbsoluteTolerance()          self.setAbsoluteTolerance()
# Line 1832  class HomogeneousSaddlePointProblem(obje Line 1832  class HomogeneousSaddlePointProblem(obje
1832               correction_step+=1               correction_step+=1
1833               v,p =v2, p2               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.3851  
changed lines
  Added in v.3852

  ViewVC Help
Powered by ViewVC 1.1.26