Diff of /trunk/escript/py_src/pdetools.py

trunk/esys2/escript/py_src/pdetools.py revision 148 by jgs, Tue Aug 23 01:24:31 2005 UTC trunk/escript/py_src/pdetools.py revision 877 by gross, Wed Oct 25 03:06:58 2006 UTC
# Line 1  Line 1
1  # \$Id\$  # \$Id\$
2
3  """  """
4  provides a some tools related to PDEs currently includes:  Provides some tools related to PDEs.
5
6       Projector - to project a discontinuous  Currently includes:
7        - 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
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"
22                        http://www.access.edu.au
26    __url__="http://www.iservo.edu.au/esys"
27    __version__="\$Revision\$"
28    __date__="\$Date\$"
29
"""
30
31  from escript import *  import escript
32  from linearPDEs import LinearPDE  import linearPDEs
33  import numarray  import numarray
34    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    """The Projector is a factory which projects a discontiuous function onto a continuous function on the a given domain"""    """
116      The Projector is a factory which projects a discontiuous function onto a
117      continuous function on the a given domain.
118      """
119    def __init__(self, domain, reduce = True, fast=True):    def __init__(self, domain, reduce = True, fast=True):
120      """      """
121      @brief Create a continuous function space projector for a domain.      Create a continuous function space projector for a domain.
122
123      @param domain Domain of the projection.      @param domain: Domain of the projection.
124      @param reduce Flag to reduce projection order (default is True)      @param reduce: Flag to reduce projection order (default is True)
125      @param fast Flag to use a fast method based on matrix lumping (default is true)      @param fast: Flag to use a fast method based on matrix lumping (default is true)
126      """      """
127      self.__pde = LinearPDE(domain)      self.__pde = linearPDEs.LinearPDE(domain)
128      if fast:      if fast:
129        self.__pde.setSolverMethod(LinearPDE.LUMPING)        self.__pde.setSolverMethod(linearPDEs.LinearPDE.LUMPING)
130      self.__pde.setSymmetryOn()      self.__pde.setSymmetryOn()
131      self.__pde.setReducedOrderTo(reduce)      self.__pde.setReducedOrderTo(reduce)
132      self.__pde.setValue(D = 1.)      self.__pde.setValue(D = 1.)
# Line 35  class Projector: Line 137  class Projector:
137
138    def __call__(self, input_data):    def __call__(self, input_data):
139      """      """
140      @brief projects input_data onto a continuous function      Projects input_data onto a continuous function
141
142      @param input_data  The input_data to be projected.      @param input_data: The input_data to be projected.
143      """      """
144      out=Data(0.,input_data.getShape(),what=ContinuousFunction(self.__pde.getDomain()))      out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())
145      if input_data.getRank()==0:      if input_data.getRank()==0:
146          self.__pde.setValue(Y = input_data)          self.__pde.setValue(Y = input_data)
147          out=self.__pde.getSolution()          out=self.__pde.getSolution()
# Line 67  class Projector: Line 169  class Projector:
169                      out[i0,i1,i2,i3]=self.__pde.getSolution()                      out[i0,i1,i2,i3]=self.__pde.getSolution()
170      return out      return out
171
172    class NoPDE:
class Locator:
173       """       """
174         solves the following problem for u:
175
176         M{kronecker[i,j]*D[j]*u[j]=Y[i]}
177
178         with constraint
179
180         M{u[j]=r[j]}  where M{q[j]>0}
181
182         where D, Y, r and q are given functions of rank 1.
183
184          Locator provides a  access the values of data objects at a given spatial coordinate x.       In the case of scalars this takes the form
In fact, a Locator object finds the sample in the set of samples of a given function space or domain where which is closest
to the given point x.
185
186         M{D*u=Y}
187
188         with constraint
189
190         M{u=r}  where M{q>0}
191
192         where D, Y, r and q are given scalar functions.
193
194         The constraint is overwriting any other condition.
195
196         @note: This class is similar to the L{linearPDEs.LinearPDE} class with A=B=C=X=0 but has the intention
197                that all input parameter are given in L{Solution} or L{ReducedSolution}. The whole
198                thing is a bit strange and I blame Robert.Woodcock@csiro.au for this.
199         """
200         def __init__(self,domain,D=None,Y=None,q=None,r=None):
201             """
202             initialize the problem
203
204             @param domain: domain of the PDE.
205             @type domain: L{Domain}
206             @param D: coefficient of the solution.
207             @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}
208             @param Y: right hand side
209             @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}
210             @param q: location of constraints
211             @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}
212             @param r: value of solution at locations of constraints
213             @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}
214             """
215             self.__domain=domain
216             self.__D=D
217             self.__Y=Y
218             self.__q=q
219             self.__r=r
220             self.__u=None
221             self.__function_space=escript.Solution(self.__domain)
222         def setReducedOn(self):
223             """
224             sets the L{FunctionSpace} of the solution to L{ReducedSolution}
225             """
226             self.__function_space=escript.ReducedSolution(self.__domain)
227             self.__u=None
228
229         def setReducedOff(self):
230             """
231             sets the L{FunctionSpace} of the solution to L{Solution}
232             """
233             self.__function_space=escript.Solution(self.__domain)
234             self.__u=None
235
236         def setValue(self,D=None,Y=None,q=None,r=None):
237             """
238             assigns values to the parameters.
239
240             @param D: coefficient of the solution.
241             @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}
242             @param Y: right hand side
243             @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}
244             @param q: location of constraints
245             @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}
246             @param r: value of solution at locations of constraints
247             @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}
248             """
249             if not D==None:
250                self.__D=D
251                self.__u=None
252             if not Y==None:
253                self.__Y=Y
254                self.__u=None
255             if not q==None:
256                self.__q=q
257                self.__u=None
258             if not r==None:
259                self.__r=r
260                self.__u=None
261
262         def getSolution(self):
263             """
264             returns the solution
265
266             @return: the solution of the problem
267             @rtype: L{Data} object in the L{FunctionSpace} L{Solution} or L{ReducedSolution}.
268             """
269             if self.__u==None:
270                if self.__D==None:
271                   raise ValueError,"coefficient D is undefined"
272                D=escript.Data(self.__D,self.__function_space)
273                if D.getRank()>1:
274                   raise ValueError,"coefficient D must have rank 0 or 1"
275                if self.__Y==None:
276                   self.__u=escript.Data(0.,D.getShape(),self.__function_space)
277                else:
278                   self.__u=util.quotient(self.__Y,D)
279                if not self.__q==None:
280                    q=util.wherePositive(escript.Data(self.__q,self.__function_space))
281                    self.__u*=(1.-q)
282                    if not self.__r==None: self.__u+=q*self.__r
283             return self.__u
284
285    class Locator:
286         """
287         Locator provides access to the values of data objects at a given
288         spatial coordinate x.
289
290         In fact, a Locator object finds the sample in the set of samples of a
291         given function space or domain where which is closest to the given
292         point x.
293       """       """
294
295       def __init__(self,where,x=numarray.zeros((3,))):       def __init__(self,where,x=numarray.zeros((3,))):
296         """initializes a Locator to access values in Data objects on the Doamin or FunctionSpace where for the sample point which         """
297            closest to the given point x"""         Initializes a Locator to access values in Data objects on the Doamin
298         if isinstance(where,FunctionSpace):         or FunctionSpace where for the sample point which
299           closest to the given point x.
300           """
301           if isinstance(where,escript.FunctionSpace):
302            self.__function_space=where            self.__function_space=where
303         else:         else:
304            self.__function_space=ContinuousFunction(where)            self.__function_space=escript.ContinuousFunction(where)
305         self.__id=length(x[:self.__function_space.getDim()]-self.__function_space.getX()).mindp()         self.__id=util.length(self.__function_space.getX()-x[:self.__function_space.getDim()]).mindp()
306
307       def __str__(self):       def __str__(self):
308         """returns the coordinates of the Locator as a string"""         """
309           Returns the coordinates of the Locator as a string.
310           """
311         return "<Locator %s>"%str(self.getX())         return "<Locator %s>"%str(self.getX())
312
313       def getFunctionSpace(self):       def getFunctionSpace(self):
314          """returns the function space of the Locator"""          """
315        Returns the function space of the Locator.
316        """
317          return self.__function_space          return self.__function_space
318
319       def getId(self):       def getId(self):
320          """returns the identifier of the location"""          """
321        Returns the identifier of the location.
322        """
323          return self.__id          return self.__id
324
325       def getX(self):       def getX(self):
326          """returns the exact coordinates of the Locator"""          """
327        Returns the exact coordinates of the Locator.
328        """
329          return self(self.getFunctionSpace().getX())          return self(self.getFunctionSpace().getX())
330
331       def __call__(self,data):       def __call__(self,data):
332          """returns the value of data at the Locator of a Data object otherwise the object is returned."""          """
333        Returns the value of data at the Locator of a Data object otherwise
334        the object is returned.
335        """
336          return self.getValue(data)          return self.getValue(data)
337
338       def getValue(self,data):       def getValue(self,data):
339          """returns the value of data at the Locator if data is a Data object otherwise the object is returned."""          """
340          if isinstance(data,Data):      Returns the value of data at the Locator if data is a Data object
341        otherwise the object is returned.
342        """
343            if isinstance(data,escript.Data):
344             if data.getFunctionSpace()==self.getFunctionSpace():             if data.getFunctionSpace()==self.getFunctionSpace():
345               out=data.convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])               #out=data.convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])
346                 out=data.convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1],self.getId()[2])
347             else:             else:
348               out=data.interpolate(self.getFunctionSpace()).convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])               #out=data.interpolate(self.getFunctionSpace()).convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])
349                 out=data.interpolate(self.getFunctionSpace()).convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1],self.getId()[2])
350             if data.getRank()==0:             if data.getRank()==0:
351                return out[0]                return out[0]
352             else:             else:
353                return out                return out
354          else:          else:
355             return data             return data
356
358       """
359       This implements a solver for a saddlepoint problem
360
361       M{f(u,p)=0}
362       M{g(u)=0}
363
364       for u and p. The problem is solved with an inexact Uszawa scheme for p:
365
366       M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})
367       M{Q_g (p^{k+1}-p^{k}) =   g(u^{k+1})}
368
369       where Q_f is an approximation of the Jacobiean A_f of f with respect to u  and Q_f is an approximation of
370       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'
371       Q_g can be difficult, non-linear conjugate gradient method is applied to solve for p, so Q_g plays
372       in fact the role of a preconditioner.
373       """
374       def __init__(self,verbose=False,*args):
375           """
376           initializes the problem
377
378           @parm verbose: switches on the printing out some information
379           @type verbose: C{bool}
380           @note: this method may be overwritten by a particular saddle point problem
381           """
382           self.__verbose=verbose
383           self.relaxation=1.
384
385       def trace(self,text):
386           """
387           prints text if verbose has been set
388
389           @parm text: a text message
390           @type text: C{str}
391           """
392           if self.__verbose: print "%s: %s"%(str(self),text)
393
394       def solve_f(self,u,p,tol=1.e-8):
395           """
396           solves
397
398           A_f du = f(u,p)
399
400           with tolerance C{tol} and return du. A_f is Jacobiean of f with respect to u.
401
402           @param u: current approximation of u
403           @type u: L{escript.Data}
404           @param p: current approximation of p
405           @type p: L{escript.Data}
406           @param tol: tolerance expected for du
407           @type tol: C{float}
408           @return: increment du
409           @rtype: L{escript.Data}
410           @note: this method has to be overwritten by a particular saddle point problem
411           """
412           pass
413
414       def solve_g(self,u,tol=1.e-8):
415           """
416           solves
417
418           Q_g dp = g(u)
419
420           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.
421
422           @param u: current approximation of u
423           @type u: L{escript.Data}
424           @param tol: tolerance expected for dp
425           @type tol: C{float}
426           @return: increment dp
427           @rtype: L{escript.Data}
428           @note: this method has to be overwritten by a particular saddle point problem
429           """
430           pass
431
432       def inner(self,p0,p1):
433           """
434           inner product of p0 and p1 approximating p. Typically this returns integrate(p0*p1)
435           @return: inner product of p0 and p1
436           @rtype: C{float}
437           """
438           pass
439
440       subiter_max=3
441       def solve(self,u0,p0,tolerance=1.e-6,tolerance_u=None,iter_max=100,accepted_reduction=0.995,relaxation=None):
442            """
443            runs the solver
444
445            @param u0: initial guess for C{u}
446            @type u0: L{esys.escript.Data}
447            @param p0: initial guess for C{p}
448            @type p0: L{esys.escript.Data}
449            @param tolerance: tolerance for relative error in C{u} and C{p}
450            @type tolerance: positive C{float}
451            @param tolerance_u: tolerance for relative error in C{u} if different from C{tolerance}
452            @type tolerance_u: positive C{float}
453            @param iter_max: maximum number of iteration steps.
454            @type iter_max: C{int}
455            @param accepted_reduction: if the norm  g cannot be reduced by C{accepted_reduction} backtracking to adapt the
456                                       relaxation factor. If C{accepted_reduction=None} no backtracking is used.
457            @type accepted_reduction: positive C{float} or C{None}
458            @param relaxation: initial relaxation factor. If C{relaxation==None}, the last relaxation factor is used.
459            @type relaxation: C{float} or C{None}
460            """
461            tol=1.e-2
462            if tolerance_u==None: tolerance_u=tolerance
463            if not relaxation==None: self.relaxation=relaxation
464            if accepted_reduction ==None:
465                  angle_limit=0.
466            elif accepted_reduction>=1.:
467                  angle_limit=0.
468            else:
469                  angle_limit=util.sqrt(1-accepted_reduction**2)
470            self.iter=0
471            u=u0
472            p=p0
473            #
474            #   initialize things:
475            #
476            converged=False
477            #
478            #  start loop:
479            #
480            #  initial search direction is g
481            #
482            while not converged :
483                if self.iter>iter_max:
484                    raise ArithmeticError("no convergence after %s steps."%self.iter)
485                f_new=self.solve_f(u,p,tol)
486                norm_f_new = util.Lsup(f_new)
487                u_new=u-f_new
488                g_new=self.solve_g(u_new,tol)
489                self.iter+=1
490                norm_g_new = util.sqrt(self.inner(g_new,g_new))
491                if norm_f_new==0. and norm_g_new==0.: return u, p
492                if self.iter>1 and not accepted_reduction==None:
493                   #
494                   #   did we manage to reduce the norm of G? I
495                   #   if not we start a backtracking procedure
496                   #
497                   # print "new/old norm = ",norm_g_new, norm_g, norm_g_new/norm_g
498                   if norm_g_new > accepted_reduction * norm_g:
499                      sub_iter=0
500                      s=self.relaxation
501                      d=g
502                      g_last=g
503                      self.trace("    start substepping: f = %s, g = %s, relaxation = %s."%(norm_f_new, norm_g_new, s))
504                      while sub_iter < self.subiter_max and  norm_g_new > accepted_reduction * norm_g:
505                         dg= g_new-g_last
506                         norm_dg=abs(util.sqrt(self.inner(dg,dg))/self.relaxation)
510                         if abs(rad) < norm_dg*norm_g_new * angle_limit:
511                             if sub_iter>0: self.trace("    no further improvements expected from backtracking.")
512                             break
513                         r=self.relaxation
515                         s+=self.relaxation
516                         #####
517                         # a=g_new+self.relaxation*dg/r
518                         # print "predicted new norm = ",util.sqrt(self.inner(a,a)),util.sqrt(self.inner(g_new,g_new)), self.relaxation
519                         #####
520                         g_last=g_new
521                         p+=self.relaxation*d
522                         f_new=self.solve_f(u,p,tol)
523                         u_new=u-f_new
524                         g_new=self.solve_g(u_new,tol)
525                         self.iter+=1
526                         norm_f_new = util.Lsup(f_new)
527                         norm_g_new = util.sqrt(self.inner(g_new,g_new))
528                         # print "   ",sub_iter," new g norm",norm_g_new
529                         self.trace("    %s th sub-step: f = %s, g = %s, relaxation = %s."%(sub_iter, norm_f_new, norm_g_new, s))
530                         #
531                         #   can we expect reduction of g?
532                         #
533                         # u_last=u_new
534                         sub_iter+=1
535                      self.relaxation=s
536                #
537                #  check for convergence:
538                #
539                norm_u_new = util.Lsup(u_new)
540                p_new=p+self.relaxation*g_new
541                norm_p_new = util.sqrt(self.inner(p_new,p_new))
542                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))
543
544                if self.iter>1:
545                   dg2=g_new-g
546                   df2=f_new-f
547                   norm_dg2=util.sqrt(self.inner(dg2,dg2))
548                   norm_df2=util.Lsup(df2)
549                   # print norm_g_new, norm_g, norm_dg, norm_p, tolerance
550                   tol_eq_g=tolerance*norm_dg2/(norm_g*abs(self.relaxation))*norm_p_new
551                   tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new
552                   if norm_g_new <= tol_eq_g and norm_f_new <= tol_eq_f:
553                       converged=True
554                       break
555                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
556            self.trace("convergence after %s steps."%self.iter)
557            return u,p
558    #   def solve(self,u0,p0,tolerance=1.e-6,iter_max=10,self.relaxation=1.):
559    #      tol=1.e-2
560    #      iter=0
561    #      converged=False
562    #      u=u0*1.
563    #      p=p0*1.
564    #      while not converged and iter<iter_max:
565    #          du=self.solve_f(u,p,tol)
566    #          u-=du
567    #          norm_du=util.Lsup(du)
568    #          norm_u=util.Lsup(u)
569    #
570    #          dp=self.relaxation*self.solve_g(u,tol)
571    #          p+=dp
572    #          norm_dp=util.sqrt(self.inner(dp,dp))
573    #          norm_p=util.sqrt(self.inner(p,p))
574    #          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)
575    #          iter+=1
576    #
577    #          converged = (norm_du <= tolerance*norm_u) and  (norm_dp <= tolerance*norm_p)
578    #      if converged:
579    #          print "convergence after %s steps."%iter
580    #      else:
581    #          raise ArithmeticError("no convergence after %s steps."%iter)
582    #
583    #      return u,p
584
585    # vim: expandtab shiftwidth=4:

Legend:
 Removed from v.148 changed lines Added in v.877