/[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

trunk/esys2/escript/py_src/pdetools.py revision 149 by jgs, Thu Sep 1 03:31:39 2005 UTC trunk/escript/py_src/pdetools.py revision 1107 by gross, Thu Apr 19 02:14:18 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      if input_data.getRank()==0:      if input_data.getRank()==0:
143          self.__pde.setValue(Y = input_data)          self.__pde.setValue(Y = input_data)
144          out=self.__pde.getSolution()          out=self.__pde.getSolution()
# Line 70  class Projector: Line 166  class Projector:
166                      out[i0,i1,i2,i3]=self.__pde.getSolution()                      out[i0,i1,i2,i3]=self.__pde.getSolution()
167      return out      return out
168    
169    class NoPDE:
170         """
171         solves the following problem for u:
172    
173         M{kronecker[i,j]*D[j]*u[j]=Y[i]}
174    
175         with constraint
176    
177         M{u[j]=r[j]}  where M{q[j]>0}
178    
179         where D, Y, r and q are given functions of rank 1.
180    
181         In the case of scalars this takes the form
182    
183         M{D*u=Y}
184    
185         with constraint
186    
187         M{u=r}  where M{q>0}
188    
189         where D, Y, r and q are given scalar functions.
190    
191         The constraint is overwriting any other condition.
192    
193         @note: This class is similar to the L{linearPDEs.LinearPDE} class with A=B=C=X=0 but has the intention
194                that all input parameter are given in L{Solution} or L{ReducedSolution}. The whole
195                thing is a bit strange and I blame Robert.Woodcock@csiro.au for this.
196         """
197         def __init__(self,domain,D=None,Y=None,q=None,r=None):
198             """
199             initialize the problem
200    
201             @param domain: domain of the PDE.
202             @type domain: L{Domain}
203             @param D: coefficient of the solution.
204             @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}
205             @param Y: right hand side
206             @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}
207             @param q: location of constraints
208             @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}
209             @param r: value of solution at locations of constraints
210             @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}
211             """
212             self.__domain=domain
213             self.__D=D
214             self.__Y=Y
215             self.__q=q
216             self.__r=r
217             self.__u=None
218             self.__function_space=escript.Solution(self.__domain)
219         def setReducedOn(self):
220             """
221             sets the L{FunctionSpace} of the solution to L{ReducedSolution}
222             """
223             self.__function_space=escript.ReducedSolution(self.__domain)
224             self.__u=None
225    
226         def setReducedOff(self):
227             """
228             sets the L{FunctionSpace} of the solution to L{Solution}
229             """
230             self.__function_space=escript.Solution(self.__domain)
231             self.__u=None
232            
233         def setValue(self,D=None,Y=None,q=None,r=None):
234             """
235             assigns values to the parameters.
236    
237             @param D: coefficient of the solution.
238             @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}
239             @param Y: right hand side
240             @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}
241             @param q: location of constraints
242             @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}
243             @param r: value of solution at locations of constraints
244             @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}
245             """
246             if not D==None:
247                self.__D=D
248                self.__u=None
249             if not Y==None:
250                self.__Y=Y
251                self.__u=None
252             if not q==None:
253                self.__q=q
254                self.__u=None
255             if not r==None:
256                self.__r=r
257                self.__u=None
258    
259         def getSolution(self):
260             """
261             returns the solution
262            
263             @return: the solution of the problem
264             @rtype: L{Data} object in the L{FunctionSpace} L{Solution} or L{ReducedSolution}.
265             """
266             if self.__u==None:
267                if self.__D==None:
268                   raise ValueError,"coefficient D is undefined"
269                D=escript.Data(self.__D,self.__function_space)
270                if D.getRank()>1:
271                   raise ValueError,"coefficient D must have rank 0 or 1"
272                if self.__Y==None:
273                   self.__u=escript.Data(0.,D.getShape(),self.__function_space)
274                else:
275                   self.__u=util.quotient(self.__Y,D)
276                if not self.__q==None:
277                    q=util.wherePositive(escript.Data(self.__q,self.__function_space))
278                    self.__u*=(1.-q)
279                    if not self.__r==None: self.__u+=q*self.__r
280             return self.__u
281                
282  class Locator:  class Locator:
283       """       """
284       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 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(x[:self.__function_space.getDim()]-self.__function_space.getX()).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 105  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 131  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
373             else:             else:
374               out=data.interpolate(self.getFunctionSpace()).convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])               dat=data.interpolate(self.getFunctionSpace())
375             if data.getRank()==0:             id=self.getId()
376                return out[0]             r=data.getRank()
377               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.149  
changed lines
  Added in v.1107

  ViewVC Help
Powered by ViewVC 1.1.26