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

Legend:
Removed from v.790  
changed lines
  Added in v.1107

  ViewVC Help
Powered by ViewVC 1.1.26