/[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 351 by gross, Tue Dec 13 09:12:15 2005 UTC revision 877 by gross, Wed Oct 25 03:06:58 2006 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
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
# Line 18  class TimeIntegrationManager: Line 37  class TimeIntegrationManager:
37    """    """
38    a simple mechanism to manage time dependend values.    a simple mechanism to manage time dependend values.
39    
40    typical usage is:    typical usage is::
41    
42    dt=0.1 # time increment       dt=0.1 # time increment
43    tm=TimeIntegrationManager(inital_value,p=1)       tm=TimeIntegrationManager(inital_value,p=1)
44    while t<1.       while t<1.
45        v_guess=tm.extrapolate(dt) # extrapolate to t+dt           v_guess=tm.extrapolate(dt) # extrapolate to t+dt
46        v=...           v=...
47        tm.checkin(dt,v)           tm.checkin(dt,v)
48        t+=dt           t+=dt
49    
50    @remark: currently only p=1 is supported.    @note: currently only p=1 is supported.
51    """    """
52    def __init__(self,*inital_values,**kwargs):    def __init__(self,*inital_values,**kwargs):
53       """       """
# Line 49  class TimeIntegrationManager: Line 68  class TimeIntegrationManager:
68    
69    def getTime(self):    def getTime(self):
70        return self.__t        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):    def checkin(self,dt,*values):
79        """        """
# Line 84  class TimeIntegrationManager: Line 109  class TimeIntegrationManager:
109           return out[0]           return out[0]
110        else:        else:
111           return out           return out
112    
113    
114  class Projector:  class Projector:
115    """    """
# Line 115  class Projector: Line 141  class Projector:
141    
142      @param input_data: The input_data to be projected.      @param input_data: The input_data to be projected.
143      """      """
144      out=escript.Data(0.,input_data.getShape(),what=escript.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 143  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:
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         In the case of scalars this takes the form
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:  class Locator:
286       """       """
287       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 164  class Locator: Line 302  class Locator:
302            self.__function_space=where            self.__function_space=where
303         else:         else:
304            self.__function_space=escript.ContinuousFunction(where)            self.__function_space=escript.ContinuousFunction(where)
305         self.__id=util.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         """         """
# Line 204  class Locator: Line 342  class Locator:
342      """      """
343          if isinstance(data,escript.Data):          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:
# Line 214  class Locator: Line 354  class Locator:
354          else:          else:
355             return data             return data
356    
357    class SaddlePointProblem(object):
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)
507                         rad=self.inner(g_new,dg)/self.relaxation
508                         # print "   ",sub_iter,": rad, norm_dg:",abs(rad), norm_dg*norm_g_new * angle_limit
509                         # print "   ",sub_iter,": rad, norm_dg:",rad, norm_dg, norm_g_new, norm_g
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
514                         self.relaxation= - rad/norm_dg**2
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:  # vim: expandtab shiftwidth=4:

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

  ViewVC Help
Powered by ViewVC 1.1.26