/[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 146 by jgs, Fri Jul 29 01:44:29 2005 UTC trunk/escript/py_src/pdetools.py revision 782 by bcumming, Tue Jul 18 00:47:47 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    
11    @var __author__: name of author
12    @var __copyright__: copyrights
13    @var __license__: licence agreement
14    @var __url__: url entry point on documentation
15    @var __version__: version
16    @var __date__: date of the version
17    """
18    
19    __author__="Lutz Gross, l.gross@uq.edu.au"
20    __copyright__="""  Copyright (c) 2006 by ACcESS MNRF
21                        http://www.access.edu.au
22                    Primary Business: Queensland, Australia"""
23    __license__="""Licensed under the Open Software License version 3.0
24                 http://www.opensource.org/licenses/osl-3.0.php"""
25    __url__="http://www.iservo.edu.au/esys"
26    __version__="$Revision$"
27    __date__="$Date$"
28    
 """  
29    
30  from escript import *  import escript
31  from linearPDEs import LinearPDE  import linearPDEs
32  import numarray  import numarray
33    import util
34    
35    class TimeIntegrationManager:
36      """
37      a simple mechanism to manage time dependend values.
38    
39      typical usage is::
40    
41         dt=0.1 # time increment
42         tm=TimeIntegrationManager(inital_value,p=1)
43         while t<1.
44             v_guess=tm.extrapolate(dt) # extrapolate to t+dt
45             v=...
46             tm.checkin(dt,v)
47             t+=dt
48    
49      @note: currently only p=1 is supported.
50      """
51      def __init__(self,*inital_values,**kwargs):
52         """
53         sets up the value manager where inital_value is the initial value and p is order used for extrapolation
54         """
55         if kwargs.has_key("p"):
56                self.__p=kwargs["p"]
57         else:
58                self.__p=1
59         if kwargs.has_key("time"):
60                self.__t=kwargs["time"]
61         else:
62                self.__t=0.
63         self.__v_mem=[inital_values]
64         self.__order=0
65         self.__dt_mem=[]
66         self.__num_val=len(inital_values)
67    
68      def getTime(self):
69          return self.__t
70      def getValue(self):
71          out=self.__v_mem[0]
72          if len(out)==1:
73              return out[0]
74          else:
75              return out
76    
77      def checkin(self,dt,*values):
78          """
79          adds new values to the manager. the p+1 last value get lost
80          """
81          o=min(self.__order+1,self.__p)
82          self.__order=min(self.__order+1,self.__p)
83          v_mem_new=[values]
84          dt_mem_new=[dt]
85          for i in range(o-1):
86             v_mem_new.append(self.__v_mem[i])
87             dt_mem_new.append(self.__dt_mem[i])
88          v_mem_new.append(self.__v_mem[o-1])
89          self.__order=o
90          self.__v_mem=v_mem_new
91          self.__dt_mem=dt_mem_new
92          self.__t+=dt
93    
94      def extrapolate(self,dt):
95          """
96          extrapolates to dt forward in time.
97          """
98          if self.__order==0:
99             out=self.__v_mem[0]
100          else:
101            out=[]
102            for i in range(self.__num_val):
103               out.append((1.+dt/self.__dt_mem[0])*self.__v_mem[0][i]-dt/self.__dt_mem[0]*self.__v_mem[1][i])
104    
105          if len(out)==0:
106             return None
107          elif len(out)==1:
108             return out[0]
109          else:
110             return out
111    
112    
113  class Projector:  class Projector:
114    """The Projector is a factory which projects a discontiuous function onto a continuous function on the a given domain"""    """
115      The Projector is a factory which projects a discontiuous function onto a
116      continuous function on the a given domain.
117      """
118    def __init__(self, domain, reduce = True, fast=True):    def __init__(self, domain, reduce = True, fast=True):
119      """      """
120      @brief Create a continuous function space projector for a domain.      Create a continuous function space projector for a domain.
121    
122      @param domain Domain of the projection.      @param domain: Domain of the projection.
123      @param reduce Flag to reduce projection order (default is True)      @param reduce: Flag to reduce projection order (default is True)
124      @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)
125      """      """
126      self.__pde = LinearPDE(domain)      self.__pde = linearPDEs.LinearPDE(domain)
127      self.__pde.setLumping(fast)      if fast:
128          self.__pde.setSolverMethod(linearPDEs.LinearPDE.LUMPING)
129      self.__pde.setSymmetryOn()      self.__pde.setSymmetryOn()
130      self.__pde.setReducedOrderTo(reduce)      self.__pde.setReducedOrderTo(reduce)
131      self.__pde.setValue(D = 1.)      self.__pde.setValue(D = 1.)
# Line 34  class Projector: Line 136  class Projector:
136    
137    def __call__(self, input_data):    def __call__(self, input_data):
138      """      """
139      @brief projects input_data onto a continuous function      Projects input_data onto a continuous function
140    
141      @param input_data  The input_data to be projected.      @param input_data: The input_data to be projected.
142      """      """
143      out=Data(0.,input_data.getShape(),what=ContinuousFunction(self.__pde.getDomain()))      out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())
144      if input_data.getRank()==0:      if input_data.getRank()==0:
145          self.__pde.setValue(Y = input_data)          self.__pde.setValue(Y = input_data)
146          out=self.__pde.getSolution()          out=self.__pde.getSolution()
# Line 66  class Projector: Line 168  class Projector:
168                      out[i0,i1,i2,i3]=self.__pde.getSolution()                      out[i0,i1,i2,i3]=self.__pde.getSolution()
169      return out      return out
170    
171    class NoPDE:
172  class Location:       """
173       """Location provides a factory to access the values of data objects at a given spatial coordinate x.       solves the following problem for u:
174          In fact, a Location object finds the sample in the set of samples of a given function space or domain which is closest  
175          to the given point x. """       M{kronecker[i,j]*D[j]*u[j]=Y[i]}
176    
177       def __init__(self,x,where):       with constraint
178         """initializes a Location to access values in Data objects on the Doamin or FunctionSpace where for the sample point which  
179            closest to the given point x"""       M{u[j]=r[j]}  where M{q[j]>0}
180         if isinstance(where,FunctionSpace):  
181            self.__where=where       where D, Y, r and q are given functions of rank 1.
182    
183         In the case of scalars this takes the form
184    
185         M{D*u=Y}
186    
187         with constraint
188    
189         M{u=r}  where M{q>0}
190    
191         where D, Y, r and q are given scalar functions.
192    
193         The constraint is overwriting any other condition.
194    
195         @note: This class is similar to the L{linearPDEs.LinearPDE} class with A=B=C=X=0 but has the intention
196                that all input parameter are given in L{Solution} or L{ReducedSolution}. The whole
197                thing is a bit strange and I blame Robert.Woodcock@csiro.au for this.
198         """
199         def __init__(self,domain,D=None,Y=None,q=None,r=None):
200             """
201             initialize the problem
202    
203             @param domain: domain of the PDE.
204             @type domain: L{Domain}
205             @param D: coefficient of the solution.
206             @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}
207             @param Y: right hand side
208             @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}
209             @param q: location of constraints
210             @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}
211             @param r: value of solution at locations of constraints
212             @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}
213             """
214             self.__domain=domain
215             self.__D=D
216             self.__Y=Y
217             self.__q=q
218             self.__r=r
219             self.__u=None
220             self.__function_space=escript.Solution(self.__domain)
221         def setReducedOn(self):
222             """
223             sets the L{FunctionSpace} of the solution to L{ReducedSolution}
224             """
225             self.__function_space=escript.ReducedSolution(self.__domain)
226             self.__u=None
227    
228         def setReducedOff(self):
229             """
230             sets the L{FunctionSpace} of the solution to L{Solution}
231             """
232             self.__function_space=escript.Solution(self.__domain)
233             self.__u=None
234            
235         def setValue(self,D=None,Y=None,q=None,r=None):
236             """
237             assigns values to the parameters.
238    
239             @param D: coefficient of the solution.
240             @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}
241             @param Y: right hand side
242             @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}
243             @param q: location of constraints
244             @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}
245             @param r: value of solution at locations of constraints
246             @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}
247             """
248             if not D==None:
249                self.__D=D
250                self.__u=None
251             if not Y==None:
252                self.__Y=Y
253                self.__u=None
254             if not q==None:
255                self.__q=q
256                self.__u=None
257             if not r==None:
258                self.__r=r
259                self.__u=None
260    
261         def getSolution(self):
262             """
263             returns the solution
264            
265             @return: the solution of the problem
266             @rtype: L{Data} object in the L{FunctionSpace} L{Solution} or L{ReducedSolution}.
267             """
268             if self.__u==None:
269                if self.__D==None:
270                   raise ValueError,"coefficient D is undefined"
271                D=escript.Data(self.__D,self.__function_space)
272                if D.getRank()>1:
273                   raise ValueError,"coefficient D must have rank 0 or 1"
274                if self.__Y==None:
275                   self.__u=escript.Data(0.,D.getShape(),self.__function_space)
276                else:
277                   self.__u=util.quotient(self.__Y,D)
278                if not self.__q==None:
279                    q=util.wherePositive(escript.Data(self.__q,self.__function_space))
280                    self.__u*=(1.-q)
281                    if not self.__r==None: self.__u+=q*self.__r
282             return self.__u
283                
284    class Locator:
285         """
286         Locator provides access to the values of data objects at a given
287         spatial coordinate x.  
288        
289         In fact, a Locator object finds the sample in the set of samples of a
290         given function space or domain where which is closest to the given
291         point x.
292         """
293    
294         def __init__(self,where,x=numarray.zeros((3,))):
295           """
296           Initializes a Locator to access values in Data objects on the Doamin
297           or FunctionSpace where for the sample point which
298           closest to the given point x.
299           """
300           if isinstance(where,escript.FunctionSpace):
301              self.__function_space=where
302         else:         else:
303            self.__where=ContinuousFunction(where)            self.__function_space=escript.ContinuousFunction(where)
304         self.__id=length(x-self.__where.getX()).minarg()         self.__id=util.length(self.__function_space.getX()-x[:self.__function_space.getDim()]).mindp()
   
      def getValue(self,data):  
         """returns the value of data at the Location at a numarray object"""  
         if isinstance(data,Data):  
            return data  
         else:  
            if not data.getFunctionSpace()==self.getFunctionSpace():  
              raise ValueError,"function space of data obejct does not match function space of Location"  
            else:  
              return data.getValue(self.getId())  
      def getX(self):  
         """returns the exact coordinates of the Location"""  
         return self.getValue(self.getFunctionSpace().getX())  
305    
306       def getId(self):       def __str__(self):
307          """returns the identifier of the location"""         """
308          return self.__id         Returns the coordinates of the Locator as a string.
309           """
310           return "<Locator %s>"%str(self.getX())
311    
312       def getFunctionSpace(self):       def getFunctionSpace(self):
313          """returns the function space of the Location"""          """
314        Returns the function space of the Locator.
315        """
316          return self.__function_space          return self.__function_space
317    
318       def __str__(self):       def getId(self):
319         """returns the coordinates of the Location as a string"""          """
320         return str(self.getX())      Returns the identifier of the location.
321        """
322  def testProjector(domain):          return self.__id
       """runs a few test of the Projector factory and returns the largest error plus a description of the test this error occured"""  
       error_max=0.  
       error_text=""  
       x=ContinuousFunction(domain).getX()  
       for f in [True,False]:  
          p=Projector(domain,reduce=False,fast=f)  
          for r in range(5):  
             text="range %s , fast=%s"%(r,f)  
             if r==0:  
                td_ref=x[0]  
             elif r==1:  
                td_ref=x  
             elif r==2:  
                td_ref=[[11.,12.],[21,22.]]*(x[0]+x[1])  
             elif r==3:  
                td_ref=[[[111.,112.],[121,122.]],[[211.,212.],[221,222.]]]*(x[0]+x[1])  
             elif r==3:  
                td_ref=[[[[1111.,1112.],[1121,1122.]],[[1211.,1212.],[1221,1222.]]],[[[2111.,2112.],[2121,2122.]],[[2211.,2212.],[2221,2222.]]]]*(x[0]+x[1])  
             td=p(td_ref.interpolate(Function(domain)))  
             er=Lsup(td-td_ref)/Lsup(td_ref)  
             print text," error = ",er  
             if er>error_max:  
                  error_max=er  
                  error_text=text  
       return error_max,error_text  
   
   
 # this should be removed later  
 if __name__=="__main__":  
     from esys.finley import Rectangle  
     txt=testProjector(Rectangle(56,61))  
     print "test Projector: ",txt  
323    
324         def getX(self):
325            """
326        Returns the exact coordinates of the Locator.
327        """
328            return self(self.getFunctionSpace().getX())
329    
330         def __call__(self,data):
331            """
332        Returns the value of data at the Locator of a Data object otherwise
333        the object is returned.
334        """
335            return self.getValue(data)
336    
337         def getValue(self,data):
338            """
339        Returns the value of data at the Locator if data is a Data object
340        otherwise the object is returned.
341        """
342            if isinstance(data,escript.Data):
343               if data.getFunctionSpace()==self.getFunctionSpace():
344                 out=data.convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])
345                 #out=data.convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1],self.getId()[2])
346               else:
347                 out=data.interpolate(self.getFunctionSpace()).convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])
348                 #out=data.interpolate(self.getFunctionSpace()).convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1],self.getId()[2])
349               if data.getRank()==0:
350                  return out[0]
351               else:
352                  return out
353            else:
354               return data
355    
356    # vim: expandtab shiftwidth=4:

Legend:
Removed from v.146  
changed lines
  Added in v.782

  ViewVC Help
Powered by ViewVC 1.1.26