/[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 873 by gross, Mon Oct 16 04:07:33 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       f(u,p)=0
362       g(u)=0
363    
364       for u and p. The problem is solved with an inexact Uszawa scheme for p:
365    
366       Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})
367       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    
384       def trace(self,text):
385           """
386           prints text if verbose has been set
387    
388           @parm text: a text message
389           @type text: C{str}
390           """
391           if self.__verbose: print "%s: %s"%(str(self),text)
392    
393       def solve_f(self,u,p,tol=1.e-8):
394           """
395           solves
396    
397           A_f du = f(u,p)
398    
399           with tolerance C{tol} and return du. A_f is Jacobiean of f with respect to u.
400    
401           @param u: current approximation of u
402           @type u: L{escript.Data}
403           @param p: current approximation of p
404           @type p: L{escript.Data}
405           @param tol: tolerance expected for du
406           @type tol: C{float}
407           @return: increment du
408           @rtype: L{escript.Data}
409           @note: this method has to be overwritten by a particular saddle point problem
410           """
411           pass
412    
413       def solve_g(self,u,tol=1.e-8):
414           """
415           solves
416    
417           Q_g dp = g(u)
418    
419           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.
420    
421           @param u: current approximation of u
422           @type u: L{escript.Data}
423           @param tol: tolerance expected for dp
424           @type tol: C{float}
425           @return: increment dp
426           @rtype: L{escript.Data}
427           @note: this method has to be overwritten by a particular saddle point problem
428           """
429           pass
430    
431       def inner(self,p0,p1):
432           """
433           inner product of p0 and p1 approximating p. Typically this returns integrate(p0*p1)
434           @return: inner product of p0 and p1
435           @rtype: C{float}
436           """
437           pass
438    
439       def solve(self,u0,p0,tolerance=1.e-6,iter_max=10,relaxation=1.):
440          tol=1.e-2
441          iter=0
442          converged=False
443          u=u0*1.
444          p=p0*1.
445          while not converged and iter<iter_max:
446              du=self.solve_f(u,p,tol)
447              u-=du
448              norm_du=util.Lsup(du)
449              norm_u=util.Lsup(u)
450            
451              dp=relaxation*self.solve_g(u,tol)
452              p+=dp
453              norm_dp=util.sqrt(self.inner(dp,dp))
454              norm_p=util.sqrt(self.inner(p,p))
455              print iter,"-th step rel. errror u,p= (%s,%s)(%s,%s)"%(norm_du/norm_u,norm_dp/norm_p,norm_u,norm_p)
456              iter+=1
457    
458              converged = (norm_du <= tolerance*norm_u) and  (norm_dp <= tolerance*norm_p)
459          if converged:
460              print "convergence after %s steps."%iter
461          else:
462              raise ArithmeticError("no convergence after %s steps."%iter)
463    
464          return u,p
465              
466  # vim: expandtab shiftwidth=4:  # vim: expandtab shiftwidth=4:

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

  ViewVC Help
Powered by ViewVC 1.1.26