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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 790 - (show annotations)
Wed Jul 26 23:12:34 2006 UTC (13 years ago) by bcumming
File MIME type: text/x-python
File size: 11391 byte(s)
changes to escript/py_src/pdetools.py and /escript/src/Data.h/.cpp to
make the Locator work in MPI. escript::Data::mindp now returns a 3 tuple,
with the MPI rank of the process on which the minimum value occurs
included. escript::Data::convertToNumArrayFromDPNo also takes the ProcNo
to perform the MPI reduction.

This had to be implemented in both the MPI and non-MPI versions to allow
the necesary changes to the Python code in pdetools.py. In the non-MPI
version ProcNo is set to 0. This works for the explicit scripts tested
thus far, however if it causes problems in your scripts contact Ben or
Lutz, or revert the three files (pdetools.py, Data.h and Data.cpp) to
the previous version.  


1 # $Id$
2
3 """
4 Provides some tools related to PDEs.
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 @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 import escript
31 import linearPDEs
32 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:
114 """
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):
119 """
120 Create a continuous function space projector for a domain.
121
122 @param domain: Domain of the projection.
123 @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)
125 """
126 self.__pde = linearPDEs.LinearPDE(domain)
127 if fast:
128 self.__pde.setSolverMethod(linearPDEs.LinearPDE.LUMPING)
129 self.__pde.setSymmetryOn()
130 self.__pde.setReducedOrderTo(reduce)
131 self.__pde.setValue(D = 1.)
132 return
133
134 def __del__(self):
135 return
136
137 def __call__(self, input_data):
138 """
139 Projects input_data onto a continuous function
140
141 @param input_data: The input_data to be projected.
142 """
143 out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())
144 if input_data.getRank()==0:
145 self.__pde.setValue(Y = input_data)
146 out=self.__pde.getSolution()
147 elif input_data.getRank()==1:
148 for i0 in range(input_data.getShape()[0]):
149 self.__pde.setValue(Y = input_data[i0])
150 out[i0]=self.__pde.getSolution()
151 elif input_data.getRank()==2:
152 for i0 in range(input_data.getShape()[0]):
153 for i1 in range(input_data.getShape()[1]):
154 self.__pde.setValue(Y = input_data[i0,i1])
155 out[i0,i1]=self.__pde.getSolution()
156 elif input_data.getRank()==3:
157 for i0 in range(input_data.getShape()[0]):
158 for i1 in range(input_data.getShape()[1]):
159 for i2 in range(input_data.getShape()[2]):
160 self.__pde.setValue(Y = input_data[i0,i1,i2])
161 out[i0,i1,i2]=self.__pde.getSolution()
162 else:
163 for i0 in range(input_data.getShape()[0]):
164 for i1 in range(input_data.getShape()[1]):
165 for i2 in range(input_data.getShape()[2]):
166 for i3 in range(input_data.getShape()[3]):
167 self.__pde.setValue(Y = input_data[i0,i1,i2,i3])
168 out[i0,i1,i2,i3]=self.__pde.getSolution()
169 return out
170
171 class NoPDE:
172 """
173 solves the following problem for u:
174
175 M{kronecker[i,j]*D[j]*u[j]=Y[i]}
176
177 with constraint
178
179 M{u[j]=r[j]} where M{q[j]>0}
180
181 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:
303 self.__function_space=escript.ContinuousFunction(where)
304 self.__id=util.length(self.__function_space.getX()-x[:self.__function_space.getDim()]).mindp()
305
306 def __str__(self):
307 """
308 Returns the coordinates of the Locator as a string.
309 """
310 return "<Locator %s>"%str(self.getX())
311
312 def getFunctionSpace(self):
313 """
314 Returns the function space of the Locator.
315 """
316 return self.__function_space
317
318 def getId(self):
319 """
320 Returns the identifier of the location.
321 """
322 return self.__id
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:

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26