/[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 867 by gross, Mon Oct 9 06:50:09 2006 UTC revision 880 by gross, Wed Oct 25 23:58:16 2006 UTC
# Line 297  class Locator: Line 297  class Locator:
297         Initializes a Locator to access values in Data objects on the Doamin         Initializes a Locator to access values in Data objects on the Doamin
298         or FunctionSpace where for the sample point which         or FunctionSpace where for the sample point which
299         closest to the given point x.         closest to the given point x.
300    
301           @param where: function space
302           @type where: L{escript.FunctionSpace}
303           @param x: coefficient of the solution.
304           @type x: L{numarray.NumArray} or C{list} of L{numarray.NumArray}
305         """         """
306         if isinstance(where,escript.FunctionSpace):         if isinstance(where,escript.FunctionSpace):
307            self.__function_space=where            self.__function_space=where
308         else:         else:
309            self.__function_space=escript.ContinuousFunction(where)            self.__function_space=escript.ContinuousFunction(where)
310         self.__id=util.length(self.__function_space.getX()-x[:self.__function_space.getDim()]).mindp()         if isinstance(x, list):
311               self.__id=[]
312               for p in x:
313                  self.__id.append(util.length(self.__function_space.getX()-p[:self.__function_space.getDim()]).mindp())
314           else:
315               self.__id=util.length(self.__function_space.getX()-x[:self.__function_space.getDim()]).mindp()
316    
317       def __str__(self):       def __str__(self):
318         """         """
319         Returns the coordinates of the Locator as a string.         Returns the coordinates of the Locator as a string.
320         """         """
321         return "<Locator %s>"%str(self.getX())         x=self.getX()
322           if instance(x,list):
323              out="["
324              first=True
325              for xx in x:
326                if not first:
327                    out+=","
328                else:
329                    first=False
330                out+=str(xx)
331              out+="]>"
332           else:
333              out=str(x)
334           return out
335    
336         def getX(self):
337            """
338        Returns the exact coordinates of the Locator.
339        """
340            return self(self.getFunctionSpace().getX())
341    
342       def getFunctionSpace(self):       def getFunctionSpace(self):
343          """          """
# Line 316  class Locator: Line 345  class Locator:
345      """      """
346          return self.__function_space          return self.__function_space
347    
348       def getId(self):       def getId(self,item=None):
349          """          """
350      Returns the identifier of the location.      Returns the identifier of the location.
351      """      """
352          return self.__id          if item == None:
353               return self.__id
354            else:
355               if isinstance(self.__id,list):
356                  return self.__id[item]
357               else:
358                  return self.__id
359    
      def getX(self):  
         """  
     Returns the exact coordinates of the Locator.  
     """  
         return self(self.getFunctionSpace().getX())  
360    
361       def __call__(self,data):       def __call__(self,data):
362          """          """
# Line 342  class Locator: Line 372  class Locator:
372      """      """
373          if isinstance(data,escript.Data):          if isinstance(data,escript.Data):
374             if data.getFunctionSpace()==self.getFunctionSpace():             if data.getFunctionSpace()==self.getFunctionSpace():
375               #out=data.convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])               dat=data
              out=data.convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1],self.getId()[2])  
376             else:             else:
377               #out=data.interpolate(self.getFunctionSpace()).convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])               dat=data.interpolate(self.getFunctionSpace())
378               out=data.interpolate(self.getFunctionSpace()).convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1],self.getId()[2])             id=self.getId()
379             if data.getRank()==0:             r=data.getRank()
380                return out[0]             if isinstance(id,list):
381                   out=[]
382                   for i in id:
383                      o=data.convertToNumArrayFromDPNo(*i)
384                      if data.getRank()==0:
385                         out.append(o[0])
386                      else:
387                         out.append(o)
388                   return out
389             else:             else:
390                return out               out=data.convertToNumArrayFromDPNo(*id)
391                 if data.getRank()==0:
392                    return out[0]
393                 else:
394                    return out
395          else:          else:
396             return data             return data
397    
# Line 358  class SaddlePointProblem(object): Line 399  class SaddlePointProblem(object):
399     """     """
400     This implements a solver for a saddlepoint problem     This implements a solver for a saddlepoint problem
401    
402     f(u,p)=0     M{f(u,p)=0}
403     g(u)=0     M{g(u)=0}
404    
405     for u and p. The problem is solved with an inexact Uszawa scheme for p:     for u and p. The problem is solved with an inexact Uszawa scheme for p:
406    
407     Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})     M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})
408     Q_g (p^{k+1}-p^{k}) =   g(u^{k+1})     M{Q_g (p^{k+1}-p^{k}) =   g(u^{k+1})}
409    
410     where Q_f is an approximation of the Jacobiean A_f of f with respect to u  and Q_f is an approximation of     where Q_f is an approximation of the Jacobiean A_f of f with respect to u  and Q_f is an approximation of
411     A_g A_f^{-1} A_g with A_g is the jacobiean of g with respect to p. As a the construction of a 'proper'     A_g A_f^{-1} A_g with A_g is the jacobiean of g with respect to p. As a the construction of a 'proper'
# Line 380  class SaddlePointProblem(object): Line 421  class SaddlePointProblem(object):
421         @note: this method may be overwritten by a particular saddle point problem         @note: this method may be overwritten by a particular saddle point problem
422         """         """
423         self.__verbose=verbose         self.__verbose=verbose
424           self.relaxation=1.
425    
426     def trace(self,text):     def trace(self,text):
427         """         """
# Line 390  class SaddlePointProblem(object): Line 432  class SaddlePointProblem(object):
432         """         """
433         if self.__verbose: print "%s: %s"%(str(self),text)         if self.__verbose: print "%s: %s"%(str(self),text)
434    
435     def solve_f(self,u,p,tol=1.e-7,*args):     def solve_f(self,u,p,tol=1.e-8):
436         """         """
437         solves         solves
438    
# Line 402  class SaddlePointProblem(object): Line 444  class SaddlePointProblem(object):
444         @type u: L{escript.Data}         @type u: L{escript.Data}
445         @param p: current approximation of p         @param p: current approximation of p
446         @type p: L{escript.Data}         @type p: L{escript.Data}
447         @param tol: tolerance for du         @param tol: tolerance expected for du
448         @type tol: C{float}         @type tol: C{float}
449         @return: increment du         @return: increment du
450         @rtype: L{escript.Data}         @rtype: L{escript.Data}
# Line 410  class SaddlePointProblem(object): Line 452  class SaddlePointProblem(object):
452         """         """
453         pass         pass
454    
455     def solve_g(self,u,*args):     def solve_g(self,u,tol=1.e-8):
456         """         """
457         solves         solves
458    
# Line 420  class SaddlePointProblem(object): Line 462  class SaddlePointProblem(object):
462    
463         @param u: current approximation of u         @param u: current approximation of u
464         @type u: L{escript.Data}         @type u: L{escript.Data}
465           @param tol: tolerance expected for dp
466           @type tol: C{float}
467         @return: increment dp         @return: increment dp
468         @rtype: L{escript.Data}         @rtype: L{escript.Data}
469         @note: this method has to be overwritten by a particular saddle point problem         @note: this method has to be overwritten by a particular saddle point problem
# Line 434  class SaddlePointProblem(object): Line 478  class SaddlePointProblem(object):
478         """         """
479         pass         pass
480    
481     def solve(self,u0,p0,tolerance=1.e-6,*args):     subiter_max=3
482         pass     def solve(self,u0,p0,tolerance=1.e-6,tolerance_u=None,iter_max=100,accepted_reduction=0.995,relaxation=None):
483            """
484            runs the solver
485    
486            @param u0: initial guess for C{u}
487            @type u0: L{esys.escript.Data}
488            @param p0: initial guess for C{p}
489            @type p0: L{esys.escript.Data}
490            @param tolerance: tolerance for relative error in C{u} and C{p}
491            @type tolerance: positive C{float}
492            @param tolerance_u: tolerance for relative error in C{u} if different from C{tolerance}
493            @type tolerance_u: positive C{float}
494            @param iter_max: maximum number of iteration steps.
495            @type iter_max: C{int}
496            @param accepted_reduction: if the norm  g cannot be reduced by C{accepted_reduction} backtracking to adapt the
497                                       relaxation factor. If C{accepted_reduction=None} no backtracking is used.
498            @type accepted_reduction: positive C{float} or C{None}
499            @param relaxation: initial relaxation factor. If C{relaxation==None}, the last relaxation factor is used.
500            @type relaxation: C{float} or C{None}
501            """
502            tol=1.e-2
503            if tolerance_u==None: tolerance_u=tolerance
504            if not relaxation==None: self.relaxation=relaxation
505            if accepted_reduction ==None:
506                  angle_limit=0.
507            elif accepted_reduction>=1.:
508                  angle_limit=0.
509            else:
510                  angle_limit=util.sqrt(1-accepted_reduction**2)
511            self.iter=0
512            u=u0
513            p=p0
514            #
515            #   initialize things:
516            #
517            converged=False
518            #
519            #  start loop:
520            #
521            #  initial search direction is g
522            #
523            while not converged :
524                if self.iter>iter_max:
525                    raise ArithmeticError("no convergence after %s steps."%self.iter)
526                f_new=self.solve_f(u,p,tol)
527                norm_f_new = util.Lsup(f_new)
528                u_new=u-f_new
529                g_new=self.solve_g(u_new,tol)
530                self.iter+=1
531                norm_g_new = util.sqrt(self.inner(g_new,g_new))
532                if norm_f_new==0. and norm_g_new==0.: return u, p
533                if self.iter>1 and not accepted_reduction==None:
534                   #
535                   #   did we manage to reduce the norm of G? I
536                   #   if not we start a backtracking procedure
537                   #
538                   # print "new/old norm = ",norm_g_new, norm_g, norm_g_new/norm_g
539                   if norm_g_new > accepted_reduction * norm_g:
540                      sub_iter=0
541                      s=self.relaxation
542                      d=g
543                      g_last=g
544                      self.trace("    start substepping: f = %s, g = %s, relaxation = %s."%(norm_f_new, norm_g_new, s))
545                      while sub_iter < self.subiter_max and  norm_g_new > accepted_reduction * norm_g:
546                         dg= g_new-g_last
547                         norm_dg=abs(util.sqrt(self.inner(dg,dg))/self.relaxation)
548                         rad=self.inner(g_new,dg)/self.relaxation
549                         # print "   ",sub_iter,": rad, norm_dg:",abs(rad), norm_dg*norm_g_new * angle_limit
550                         # print "   ",sub_iter,": rad, norm_dg:",rad, norm_dg, norm_g_new, norm_g
551                         if abs(rad) < norm_dg*norm_g_new * angle_limit:
552                             if sub_iter>0: self.trace("    no further improvements expected from backtracking.")
553                             break
554                         r=self.relaxation
555                         self.relaxation= - rad/norm_dg**2
556                         s+=self.relaxation
557                         #####
558                         # a=g_new+self.relaxation*dg/r
559                         # print "predicted new norm = ",util.sqrt(self.inner(a,a)),util.sqrt(self.inner(g_new,g_new)), self.relaxation
560                         #####
561                         g_last=g_new
562                         p+=self.relaxation*d
563                         f_new=self.solve_f(u,p,tol)
564                         u_new=u-f_new
565                         g_new=self.solve_g(u_new,tol)
566                         self.iter+=1
567                         norm_f_new = util.Lsup(f_new)
568                         norm_g_new = util.sqrt(self.inner(g_new,g_new))
569                         # print "   ",sub_iter," new g norm",norm_g_new
570                         self.trace("    %s th sub-step: f = %s, g = %s, relaxation = %s."%(sub_iter, norm_f_new, norm_g_new, s))
571                         #
572                         #   can we expect reduction of g?
573                         #
574                         # u_last=u_new
575                         sub_iter+=1
576                      self.relaxation=s
577                #
578                #  check for convergence:
579                #
580                norm_u_new = util.Lsup(u_new)
581                p_new=p+self.relaxation*g_new
582                norm_p_new = util.sqrt(self.inner(p_new,p_new))
583                self.trace("%s th step: f/u = %s, g/p = %s, relaxation = %s."%(self.iter,norm_f_new/norm_u_new, norm_g_new/norm_p_new, self.relaxation))
584    
585                if self.iter>1:
586                   dg2=g_new-g
587                   df2=f_new-f
588                   norm_dg2=util.sqrt(self.inner(dg2,dg2))
589                   norm_df2=util.Lsup(df2)
590                   # print norm_g_new, norm_g, norm_dg, norm_p, tolerance
591                   tol_eq_g=tolerance*norm_dg2/(norm_g*abs(self.relaxation))*norm_p_new
592                   tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new
593                   if norm_g_new <= tol_eq_g and norm_f_new <= tol_eq_f:
594                       converged=True
595                       break
596                f, norm_f, u, norm_u, g, norm_g, p, norm_p = f_new, norm_f_new, u_new, norm_u_new, g_new, norm_g_new, p_new, norm_p_new
597            self.trace("convergence after %s steps."%self.iter)
598            return u,p
599    #   def solve(self,u0,p0,tolerance=1.e-6,iter_max=10,self.relaxation=1.):
600    #      tol=1.e-2
601    #      iter=0
602    #      converged=False
603    #      u=u0*1.
604    #      p=p0*1.
605    #      while not converged and iter<iter_max:
606    #          du=self.solve_f(u,p,tol)
607    #          u-=du
608    #          norm_du=util.Lsup(du)
609    #          norm_u=util.Lsup(u)
610    #        
611    #          dp=self.relaxation*self.solve_g(u,tol)
612    #          p+=dp
613    #          norm_dp=util.sqrt(self.inner(dp,dp))
614    #          norm_p=util.sqrt(self.inner(p,p))
615    #          print iter,"-th step rel. errror u,p= (%s,%s) (%s,%s)(%s,%s)"%(norm_du,norm_dp,norm_du/norm_u,norm_dp/norm_p,norm_u,norm_p)
616    #          iter+=1
617    #
618    #          converged = (norm_du <= tolerance*norm_u) and  (norm_dp <= tolerance*norm_p)
619    #      if converged:
620    #          print "convergence after %s steps."%iter
621    #      else:
622    #          raise ArithmeticError("no convergence after %s steps."%iter)
623    #
624    #      return u,p
625              
626  # vim: expandtab shiftwidth=4:  # vim: expandtab shiftwidth=4:

Legend:
Removed from v.867  
changed lines
  Added in v.880

  ViewVC Help
Powered by ViewVC 1.1.26