/[escript]/temp/escript/py_src/pdetools.py
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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

Legend:
Removed from v.121  
changed lines
  Added in v.614

  ViewVC Help
Powered by ViewVC 1.1.26