/[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 155 by jgs, Wed Nov 9 02:02:19 2005 UTC revision 1122 by gross, Tue May 1 03:21:04 2007 UTC
# Line 5  Provides some tools related to PDEs. Line 5  Provides some tools related to PDEs.
5    
6  Currently includes:  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
9        - 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
13    @var __copyright__: copyrights
14    @var __license__: licence agreement
15    @var __url__: url entry point on documentation
16    @var __version__: version
17    @var __date__: date of the version
18  """  """
19    
20    __author__="Lutz Gross, l.gross@uq.edu.au"
21    __copyright__="""  Copyright (c) 2006 by ACcESS MNRF
22                        http://www.access.edu.au
23                    Primary Business: Queensland, Australia"""
24    __license__="""Licensed under the Open Software License version 3.0
25                 http://www.opensource.org/licenses/osl-3.0.php"""
26    __url__="http://www.iservo.edu.au/esys"
27    __version__="$Revision$"
28    __date__="$Date$"
29    
30    
31  import escript  import escript
32  import linearPDEs  import linearPDEs
33  import numarray  import numarray
34  import util  import util
35    
36    class TimeIntegrationManager:
37      """
38      a simple mechanism to manage time dependend values.
39    
40      typical usage is::
41    
42         dt=0.1 # time increment
43         tm=TimeIntegrationManager(inital_value,p=1)
44         while t<1.
45             v_guess=tm.extrapolate(dt) # extrapolate to t+dt
46             v=...
47             tm.checkin(dt,v)
48             t+=dt
49    
50      @note: currently only p=1 is supported.
51      """
52      def __init__(self,*inital_values,**kwargs):
53         """
54         sets up the value manager where inital_value is the initial value and p is order used for extrapolation
55         """
56         if kwargs.has_key("p"):
57                self.__p=kwargs["p"]
58         else:
59                self.__p=1
60         if kwargs.has_key("time"):
61                self.__t=kwargs["time"]
62         else:
63                self.__t=0.
64         self.__v_mem=[inital_values]
65         self.__order=0
66         self.__dt_mem=[]
67         self.__num_val=len(inital_values)
68    
69      def getTime(self):
70          return self.__t
71      def getValue(self):
72          out=self.__v_mem[0]
73          if len(out)==1:
74              return out[0]
75          else:
76              return out
77    
78      def checkin(self,dt,*values):
79          """
80          adds new values to the manager. the p+1 last value get lost
81          """
82          o=min(self.__order+1,self.__p)
83          self.__order=min(self.__order+1,self.__p)
84          v_mem_new=[values]
85          dt_mem_new=[dt]
86          for i in range(o-1):
87             v_mem_new.append(self.__v_mem[i])
88             dt_mem_new.append(self.__dt_mem[i])
89          v_mem_new.append(self.__v_mem[o-1])
90          self.__order=o
91          self.__v_mem=v_mem_new
92          self.__dt_mem=dt_mem_new
93          self.__t+=dt
94    
95      def extrapolate(self,dt):
96          """
97          extrapolates to dt forward in time.
98          """
99          if self.__order==0:
100             out=self.__v_mem[0]
101          else:
102            out=[]
103            for i in range(self.__num_val):
104               out.append((1.+dt/self.__dt_mem[0])*self.__v_mem[0][i]-dt/self.__dt_mem[0]*self.__v_mem[1][i])
105    
106          if len(out)==0:
107             return None
108          elif len(out)==1:
109             return out[0]
110          else:
111             return out
112    
113    
114  class Projector:  class Projector:
115    """    """
116    The Projector is a factory which projects a discontiuous function onto a    The Projector is a factory which projects a discontiuous function onto a
# Line 33  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
138    
139      @param input_data: The input_data to be projected.      @param input_data: The input_data to be projected.
140      """      """
141      out=escript.Data(0.,input_data.getShape(),what=escript.ContinuousFunction(self.__pde.getDomain()))      out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())
142        self.__pde.setValue(Y = escript.Data(), Y_reduced = escript.Data())
143      if input_data.getRank()==0:      if input_data.getRank()==0:
144          self.__pde.setValue(Y = input_data)          self.__pde.setValue(Y = input_data)
145          out=self.__pde.getSolution()          out=self.__pde.getSolution()
# Line 70  class Projector: Line 167  class Projector:
167                      out[i0,i1,i2,i3]=self.__pde.getSolution()                      out[i0,i1,i2,i3]=self.__pde.getSolution()
168      return out      return out
169    
170    class NoPDE:
171         """
172         solves the following problem for u:
173    
174         M{kronecker[i,j]*D[j]*u[j]=Y[i]}
175    
176         with constraint
177    
178         M{u[j]=r[j]}  where M{q[j]>0}
179    
180         where D, Y, r and q are given functions of rank 1.
181    
182         In the case of scalars this takes the form
183    
184         M{D*u=Y}
185    
186         with constraint
187    
188         M{u=r}  where M{q>0}
189    
190         where D, Y, r and q are given scalar functions.
191    
192         The constraint is overwriting any other condition.
193    
194         @note: This class is similar to the L{linearPDEs.LinearPDE} class with A=B=C=X=0 but has the intention
195                that all input parameter are given in L{Solution} or L{ReducedSolution}. The whole
196                thing is a bit strange and I blame Robert.Woodcock@csiro.au for this.
197         """
198         def __init__(self,domain,D=None,Y=None,q=None,r=None):
199             """
200             initialize the problem
201    
202             @param domain: domain of the PDE.
203             @type domain: L{Domain}
204             @param D: coefficient of the solution.
205             @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}
206             @param Y: right hand side
207             @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}
208             @param q: location of constraints
209             @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}
210             @param r: value of solution at locations of constraints
211             @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}
212             """
213             self.__domain=domain
214             self.__D=D
215             self.__Y=Y
216             self.__q=q
217             self.__r=r
218             self.__u=None
219             self.__function_space=escript.Solution(self.__domain)
220         def setReducedOn(self):
221             """
222             sets the L{FunctionSpace} of the solution to L{ReducedSolution}
223             """
224             self.__function_space=escript.ReducedSolution(self.__domain)
225             self.__u=None
226    
227         def setReducedOff(self):
228             """
229             sets the L{FunctionSpace} of the solution to L{Solution}
230             """
231             self.__function_space=escript.Solution(self.__domain)
232             self.__u=None
233            
234         def setValue(self,D=None,Y=None,q=None,r=None):
235             """
236             assigns values to the parameters.
237    
238             @param D: coefficient of the solution.
239             @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}
240             @param Y: right hand side
241             @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}
242             @param q: location of constraints
243             @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}
244             @param r: value of solution at locations of constraints
245             @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}
246             """
247             if not D==None:
248                self.__D=D
249                self.__u=None
250             if not Y==None:
251                self.__Y=Y
252                self.__u=None
253             if not q==None:
254                self.__q=q
255                self.__u=None
256             if not r==None:
257                self.__r=r
258                self.__u=None
259    
260         def getSolution(self):
261             """
262             returns the solution
263            
264             @return: the solution of the problem
265             @rtype: L{Data} object in the L{FunctionSpace} L{Solution} or L{ReducedSolution}.
266             """
267             if self.__u==None:
268                if self.__D==None:
269                   raise ValueError,"coefficient D is undefined"
270                D=escript.Data(self.__D,self.__function_space)
271                if D.getRank()>1:
272                   raise ValueError,"coefficient D must have rank 0 or 1"
273                if self.__Y==None:
274                   self.__u=escript.Data(0.,D.getShape(),self.__function_space)
275                else:
276                   self.__u=util.quotient(self.__Y,D)
277                if not self.__q==None:
278                    q=util.wherePositive(escript.Data(self.__q,self.__function_space))
279                    self.__u*=(1.-q)
280                    if not self.__r==None: self.__u+=q*self.__r
281             return self.__u
282                
283  class Locator:  class Locator:
284       """       """
285       Locator provides access to the values of data objects at a given       Locator provides access to the values of data objects at a given
# Line 86  class Locator: Line 295  class Locator:
295         Initializes a Locator to access values in Data objects on the Doamin         Initializes a Locator to access values in Data objects on the Doamin
296         or FunctionSpace where for the sample point which         or FunctionSpace where for the sample point which
297         closest to the given point x.         closest to the given point x.
298    
299           @param where: function space
300           @type where: L{escript.FunctionSpace}
301           @param x: coefficient of the solution.
302           @type x: L{numarray.NumArray} or C{list} of L{numarray.NumArray}
303         """         """
304         if isinstance(where,escript.FunctionSpace):         if isinstance(where,escript.FunctionSpace):
305            self.__function_space=where            self.__function_space=where
306         else:         else:
307            self.__function_space=escript.ContinuousFunction(where)            self.__function_space=escript.ContinuousFunction(where)
308         self.__id=util.length(x[:self.__function_space.getDim()]-self.__function_space.getX()).mindp()         if isinstance(x, list):
309               self.__id=[]
310               for p in x:
311                  self.__id.append(util.length(self.__function_space.getX()-p[:self.__function_space.getDim()]).minGlobalDataPoint())
312           else:
313               self.__id=util.length(self.__function_space.getX()-x[:self.__function_space.getDim()]).minGlobalDataPoint()
314    
315       def __str__(self):       def __str__(self):
316         """         """
317         Returns the coordinates of the Locator as a string.         Returns the coordinates of the Locator as a string.
318         """         """
319         return "<Locator %s>"%str(self.getX())         x=self.getX()
320           if instance(x,list):
321              out="["
322              first=True
323              for xx in x:
324                if not first:
325                    out+=","
326                else:
327                    first=False
328                out+=str(xx)
329              out+="]>"
330           else:
331              out=str(x)
332           return out
333    
334         def getX(self):
335            """
336        Returns the exact coordinates of the Locator.
337        """
338            return self(self.getFunctionSpace().getX())
339    
340       def getFunctionSpace(self):       def getFunctionSpace(self):
341          """          """
# Line 105  class Locator: Line 343  class Locator:
343      """      """
344          return self.__function_space          return self.__function_space
345    
346       def getId(self):       def getId(self,item=None):
347          """          """
348      Returns the identifier of the location.      Returns the identifier of the location.
349      """      """
350          return self.__id          if item == None:
351               return self.__id
352            else:
353               if isinstance(self.__id,list):
354                  return self.__id[item]
355               else:
356                  return self.__id
357    
      def getX(self):  
         """  
     Returns the exact coordinates of the Locator.  
     """  
         return self(self.getFunctionSpace().getX())  
358    
359       def __call__(self,data):       def __call__(self,data):
360          """          """
# Line 131  class Locator: Line 370  class Locator:
370      """      """
371          if isinstance(data,escript.Data):          if isinstance(data,escript.Data):
372             if data.getFunctionSpace()==self.getFunctionSpace():             if data.getFunctionSpace()==self.getFunctionSpace():
373               out=data.convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])               dat=data
374             else:             else:
375               out=data.interpolate(self.getFunctionSpace()).convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])               dat=data.interpolate(self.getFunctionSpace())
376             if data.getRank()==0:             id=self.getId()
377                return out[0]             r=data.getRank()
378               if isinstance(id,list):
379                   out=[]
380                   for i in id:
381                      o=data.getValueOfGlobalDataPoint(*i)
382                      if data.getRank()==0:
383                         out.append(o[0])
384                      else:
385                         out.append(o)
386                   return out
387             else:             else:
388                return out               out=data.getValueOfGlobalDataPoint(*id)
389                 if data.getRank()==0:
390                    return out[0]
391                 else:
392                    return out
393          else:          else:
394             return data             return data
395    
396    class SaddlePointProblem(object):
397       """
398       This implements a solver for a saddlepoint problem
399    
400       M{f(u,p)=0}
401       M{g(u)=0}
402    
403       for u and p. The problem is solved with an inexact Uszawa scheme for p:
404    
405       M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})}
406       M{Q_g (p^{k+1}-p^{k}) =   g(u^{k+1})}
407    
408       where Q_f is an approximation of the Jacobiean A_f of f with respect to u  and Q_f is an approximation of
409       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'
410       Q_g can be difficult, non-linear conjugate gradient method is applied to solve for p, so Q_g plays
411       in fact the role of a preconditioner.
412       """
413       def __init__(self,verbose=False,*args):
414           """
415           initializes the problem
416    
417           @param verbose: switches on the printing out some information
418           @type verbose: C{bool}
419           @note: this method may be overwritten by a particular saddle point problem
420           """
421           if not isinstance(verbose,bool):
422                raise TypeError("verbose needs to be of type bool.")
423           self.__verbose=verbose
424           self.relaxation=1.
425    
426       def trace(self,text):
427           """
428           prints text if verbose has been set
429    
430           @param text: a text message
431           @type text: C{str}
432           """
433           if self.__verbose: print "%s: %s"%(str(self),text)
434    
435       def solve_f(self,u,p,tol=1.e-8):
436           """
437           solves
438    
439           A_f du = f(u,p)
440    
441           with tolerance C{tol} and return du. A_f is Jacobiean of f with respect to u.
442    
443           @param u: current approximation of u
444           @type u: L{escript.Data}
445           @param p: current approximation of p
446           @type p: L{escript.Data}
447           @param tol: tolerance expected for du
448           @type tol: C{float}
449           @return: increment du
450           @rtype: L{escript.Data}
451           @note: this method has to be overwritten by a particular saddle point problem
452           """
453           pass
454    
455       def solve_g(self,u,tol=1.e-8):
456           """
457           solves
458    
459           Q_g dp = g(u)
460    
461           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.
462    
463           @param u: current approximation of u
464           @type u: L{escript.Data}
465           @param tol: tolerance expected for dp
466           @type tol: C{float}
467           @return: increment dp
468           @rtype: L{escript.Data}
469           @note: this method has to be overwritten by a particular saddle point problem
470           """
471           pass
472    
473       def inner(self,p0,p1):
474           """
475           inner product of p0 and p1 approximating p. Typically this returns integrate(p0*p1)
476           @return: inner product of p0 and p1
477           @rtype: C{float}
478           """
479           pass
480    
481       subiter_max=3
482       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                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
596            self.trace("convergence after %s steps."%self.iter)
597            return u,p
598    #   def solve(self,u0,p0,tolerance=1.e-6,iter_max=10,self.relaxation=1.):
599    #      tol=1.e-2
600    #      iter=0
601    #      converged=False
602    #      u=u0*1.
603    #      p=p0*1.
604    #      while not converged and iter<iter_max:
605    #          du=self.solve_f(u,p,tol)
606    #          u-=du
607    #          norm_du=util.Lsup(du)
608    #          norm_u=util.Lsup(u)
609    #        
610    #          dp=self.relaxation*self.solve_g(u,tol)
611    #          p+=dp
612    #          norm_dp=util.sqrt(self.inner(dp,dp))
613    #          norm_p=util.sqrt(self.inner(p,p))
614    #          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)
615    #          iter+=1
616    #
617    #          converged = (norm_du <= tolerance*norm_u) and  (norm_dp <= tolerance*norm_p)
618    #      if converged:
619    #          print "convergence after %s steps."%iter
620    #      else:
621    #          raise ArithmeticError("no convergence after %s steps."%iter)
622    #
623    #      return u,p
624              
625  # vim: expandtab shiftwidth=4:  # vim: expandtab shiftwidth=4:

Legend:
Removed from v.155  
changed lines
  Added in v.1122

  ViewVC Help
Powered by ViewVC 1.1.26