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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 351 - (hide annotations)
Tue Dec 13 09:12:15 2005 UTC (13 years, 10 months ago) by gross
File MIME type: text/x-python
File size: 6566 byte(s)
time integration manger added
1 jgs 121 # $Id$
2    
3     """
4 jgs 149 Provides some tools related to PDEs.
5 jgs 121
6 jgs 149 Currently includes:
7     - Projector - to project a discontinuous
8 gross 351 - Locator - to trace values in data objects at a certain location
9     - TimeIntegrationManager - to handel extraplotion in time
10 jgs 121 """
11    
12 jgs 149 import escript
13     import linearPDEs
14 jgs 121 import numarray
15 jgs 149 import util
16 jgs 121
17 gross 351 class TimeIntegrationManager:
18     """
19     a simple mechanism to manage time dependend values.
20    
21     typical usage is:
22    
23     dt=0.1 # time increment
24     tm=TimeIntegrationManager(inital_value,p=1)
25     while t<1.
26     v_guess=tm.extrapolate(dt) # extrapolate to t+dt
27     v=...
28     tm.checkin(dt,v)
29     t+=dt
30    
31     @remark: currently only p=1 is supported.
32     """
33     def __init__(self,*inital_values,**kwargs):
34     """
35     sets up the value manager where inital_value is the initial value and p is order used for extrapolation
36     """
37     if kwargs.has_key("p"):
38     self.__p=kwargs["p"]
39     else:
40     self.__p=1
41     if kwargs.has_key("time"):
42     self.__t=kwargs["time"]
43     else:
44     self.__t=0.
45     self.__v_mem=[inital_values]
46     self.__order=0
47     self.__dt_mem=[]
48     self.__num_val=len(inital_values)
49    
50     def getTime(self):
51     return self.__t
52    
53     def checkin(self,dt,*values):
54     """
55     adds new values to the manager. the p+1 last value get lost
56     """
57     o=min(self.__order+1,self.__p)
58     self.__order=min(self.__order+1,self.__p)
59     v_mem_new=[values]
60     dt_mem_new=[dt]
61     for i in range(o-1):
62     v_mem_new.append(self.__v_mem[i])
63     dt_mem_new.append(self.__dt_mem[i])
64     v_mem_new.append(self.__v_mem[o-1])
65     self.__order=o
66     self.__v_mem=v_mem_new
67     self.__dt_mem=dt_mem_new
68     self.__t+=dt
69    
70     def extrapolate(self,dt):
71     """
72     extrapolates to dt forward in time.
73     """
74     if self.__order==0:
75     out=self.__v_mem[0]
76     else:
77     out=[]
78     for i in range(self.__num_val):
79     out.append((1.+dt/self.__dt_mem[0])*self.__v_mem[0][i]-dt/self.__dt_mem[0]*self.__v_mem[1][i])
80    
81     if len(out)==0:
82     return None
83     elif len(out)==1:
84     return out[0]
85     else:
86     return out
87    
88 jgs 121 class Projector:
89 jgs 149 """
90     The Projector is a factory which projects a discontiuous function onto a
91     continuous function on the a given domain.
92     """
93 jgs 121 def __init__(self, domain, reduce = True, fast=True):
94     """
95 jgs 149 Create a continuous function space projector for a domain.
96 jgs 121
97 jgs 149 @param domain: Domain of the projection.
98     @param reduce: Flag to reduce projection order (default is True)
99     @param fast: Flag to use a fast method based on matrix lumping (default is true)
100 jgs 121 """
101 jgs 149 self.__pde = linearPDEs.LinearPDE(domain)
102 jgs 148 if fast:
103 jgs 149 self.__pde.setSolverMethod(linearPDEs.LinearPDE.LUMPING)
104 jgs 121 self.__pde.setSymmetryOn()
105     self.__pde.setReducedOrderTo(reduce)
106     self.__pde.setValue(D = 1.)
107     return
108    
109     def __del__(self):
110     return
111    
112     def __call__(self, input_data):
113     """
114 jgs 149 Projects input_data onto a continuous function
115 jgs 121
116 jgs 149 @param input_data: The input_data to be projected.
117 jgs 121 """
118 jgs 149 out=escript.Data(0.,input_data.getShape(),what=escript.ContinuousFunction(self.__pde.getDomain()))
119 jgs 121 if input_data.getRank()==0:
120     self.__pde.setValue(Y = input_data)
121     out=self.__pde.getSolution()
122     elif input_data.getRank()==1:
123     for i0 in range(input_data.getShape()[0]):
124     self.__pde.setValue(Y = input_data[i0])
125     out[i0]=self.__pde.getSolution()
126     elif input_data.getRank()==2:
127     for i0 in range(input_data.getShape()[0]):
128     for i1 in range(input_data.getShape()[1]):
129     self.__pde.setValue(Y = input_data[i0,i1])
130     out[i0,i1]=self.__pde.getSolution()
131     elif input_data.getRank()==3:
132     for i0 in range(input_data.getShape()[0]):
133     for i1 in range(input_data.getShape()[1]):
134     for i2 in range(input_data.getShape()[2]):
135     self.__pde.setValue(Y = input_data[i0,i1,i2])
136     out[i0,i1,i2]=self.__pde.getSolution()
137     else:
138     for i0 in range(input_data.getShape()[0]):
139     for i1 in range(input_data.getShape()[1]):
140     for i2 in range(input_data.getShape()[2]):
141     for i3 in range(input_data.getShape()[3]):
142     self.__pde.setValue(Y = input_data[i0,i1,i2,i3])
143     out[i0,i1,i2,i3]=self.__pde.getSolution()
144     return out
145    
146    
147 jgs 147 class Locator:
148     """
149 jgs 149 Locator provides access to the values of data objects at a given
150     spatial coordinate x.
151    
152     In fact, a Locator object finds the sample in the set of samples of a
153     given function space or domain where which is closest to the given
154     point x.
155 jgs 147 """
156    
157     def __init__(self,where,x=numarray.zeros((3,))):
158 jgs 149 """
159     Initializes a Locator to access values in Data objects on the Doamin
160     or FunctionSpace where for the sample point which
161     closest to the given point x.
162     """
163     if isinstance(where,escript.FunctionSpace):
164 jgs 147 self.__function_space=where
165 jgs 121 else:
166 jgs 149 self.__function_space=escript.ContinuousFunction(where)
167     self.__id=util.length(x[:self.__function_space.getDim()]-self.__function_space.getX()).mindp()
168 jgs 121
169 jgs 147 def __str__(self):
170 jgs 149 """
171     Returns the coordinates of the Locator as a string.
172     """
173 jgs 147 return "<Locator %s>"%str(self.getX())
174 jgs 121
175 jgs 147 def getFunctionSpace(self):
176 jgs 149 """
177     Returns the function space of the Locator.
178     """
179 jgs 147 return self.__function_space
180    
181 jgs 121 def getId(self):
182 jgs 149 """
183     Returns the identifier of the location.
184     """
185 jgs 121 return self.__id
186    
187 jgs 147 def getX(self):
188 jgs 149 """
189     Returns the exact coordinates of the Locator.
190     """
191 jgs 147 return self(self.getFunctionSpace().getX())
192 jgs 121
193 jgs 147 def __call__(self,data):
194 jgs 149 """
195     Returns the value of data at the Locator of a Data object otherwise
196     the object is returned.
197     """
198 jgs 147 return self.getValue(data)
199 jgs 121
200 jgs 147 def getValue(self,data):
201 jgs 149 """
202     Returns the value of data at the Locator if data is a Data object
203     otherwise the object is returned.
204     """
205     if isinstance(data,escript.Data):
206 jgs 147 if data.getFunctionSpace()==self.getFunctionSpace():
207     out=data.convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])
208     else:
209     out=data.interpolate(self.getFunctionSpace()).convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])
210     if data.getRank()==0:
211     return out[0]
212     else:
213     return out
214     else:
215     return data
216 jgs 149
217     # vim: expandtab shiftwidth=4:

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26