/[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 867 - (show annotations)
Mon Oct 9 06:50:09 2006 UTC (13 years, 1 month ago) by gross
File MIME type: text/x-python
File size: 13919 byte(s)
a class to solving sattle point problems using uszawa scheme (not functional yet)
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 - 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
32 import linearPDEs
33 import numarray
34 import util
35
36 class TimeIntegrationManager:
37 """
38 a simple mechanism to manage time dependend values.
39
40 typical usage is::
41
42 dt=0.1 # time increment
43 tm=TimeIntegrationManager(inital_value,p=1)
44 while t<1.
45 v_guess=tm.extrapolate(dt) # extrapolate to t+dt
46 v=...
47 tm.checkin(dt,v)
48 t+=dt
49
50 @note: currently only p=1 is supported.
51 """
52 def __init__(self,*inital_values,**kwargs):
53 """
54 sets up the value manager where inital_value is the initial value and p is order used for extrapolation
55 """
56 if kwargs.has_key("p"):
57 self.__p=kwargs["p"]
58 else:
59 self.__p=1
60 if kwargs.has_key("time"):
61 self.__t=kwargs["time"]
62 else:
63 self.__t=0.
64 self.__v_mem=[inital_values]
65 self.__order=0
66 self.__dt_mem=[]
67 self.__num_val=len(inital_values)
68
69 def getTime(self):
70 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):
79 """
80 adds new values to the manager. the p+1 last value get lost
81 """
82 o=min(self.__order+1,self.__p)
83 self.__order=min(self.__order+1,self.__p)
84 v_mem_new=[values]
85 dt_mem_new=[dt]
86 for i in range(o-1):
87 v_mem_new.append(self.__v_mem[i])
88 dt_mem_new.append(self.__dt_mem[i])
89 v_mem_new.append(self.__v_mem[o-1])
90 self.__order=o
91 self.__v_mem=v_mem_new
92 self.__dt_mem=dt_mem_new
93 self.__t+=dt
94
95 def extrapolate(self,dt):
96 """
97 extrapolates to dt forward in time.
98 """
99 if self.__order==0:
100 out=self.__v_mem[0]
101 else:
102 out=[]
103 for i in range(self.__num_val):
104 out.append((1.+dt/self.__dt_mem[0])*self.__v_mem[0][i]-dt/self.__dt_mem[0]*self.__v_mem[1][i])
105
106 if len(out)==0:
107 return None
108 elif len(out)==1:
109 return out[0]
110 else:
111 return out
112
113
114 class Projector:
115 """
116 The Projector is a factory which projects a discontiuous function onto a
117 continuous function on the a given domain.
118 """
119 def __init__(self, domain, reduce = True, fast=True):
120 """
121 Create a continuous function space projector for a domain.
122
123 @param domain: Domain of the projection.
124 @param reduce: Flag to reduce projection order (default is True)
125 @param fast: Flag to use a fast method based on matrix lumping (default is true)
126 """
127 self.__pde = linearPDEs.LinearPDE(domain)
128 if fast:
129 self.__pde.setSolverMethod(linearPDEs.LinearPDE.LUMPING)
130 self.__pde.setSymmetryOn()
131 self.__pde.setReducedOrderTo(reduce)
132 self.__pde.setValue(D = 1.)
133 return
134
135 def __del__(self):
136 return
137
138 def __call__(self, input_data):
139 """
140 Projects input_data onto a continuous function
141
142 @param input_data: The input_data to be projected.
143 """
144 out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())
145 if input_data.getRank()==0:
146 self.__pde.setValue(Y = input_data)
147 out=self.__pde.getSolution()
148 elif input_data.getRank()==1:
149 for i0 in range(input_data.getShape()[0]):
150 self.__pde.setValue(Y = input_data[i0])
151 out[i0]=self.__pde.getSolution()
152 elif input_data.getRank()==2:
153 for i0 in range(input_data.getShape()[0]):
154 for i1 in range(input_data.getShape()[1]):
155 self.__pde.setValue(Y = input_data[i0,i1])
156 out[i0,i1]=self.__pde.getSolution()
157 elif input_data.getRank()==3:
158 for i0 in range(input_data.getShape()[0]):
159 for i1 in range(input_data.getShape()[1]):
160 for i2 in range(input_data.getShape()[2]):
161 self.__pde.setValue(Y = input_data[i0,i1,i2])
162 out[i0,i1,i2]=self.__pde.getSolution()
163 else:
164 for i0 in range(input_data.getShape()[0]):
165 for i1 in range(input_data.getShape()[1]):
166 for i2 in range(input_data.getShape()[2]):
167 for i3 in range(input_data.getShape()[3]):
168 self.__pde.setValue(Y = input_data[i0,i1,i2,i3])
169 out[i0,i1,i2,i3]=self.__pde.getSolution()
170 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:
286 """
287 Locator provides access to the values of data objects at a given
288 spatial coordinate x.
289
290 In fact, a Locator object finds the sample in the set of samples of a
291 given function space or domain where which is closest to the given
292 point x.
293 """
294
295 def __init__(self,where,x=numarray.zeros((3,))):
296 """
297 Initializes a Locator to access values in Data objects on the Doamin
298 or FunctionSpace where for the sample point which
299 closest to the given point x.
300 """
301 if isinstance(where,escript.FunctionSpace):
302 self.__function_space=where
303 else:
304 self.__function_space=escript.ContinuousFunction(where)
305 self.__id=util.length(self.__function_space.getX()-x[:self.__function_space.getDim()]).mindp()
306
307 def __str__(self):
308 """
309 Returns the coordinates of the Locator as a string.
310 """
311 return "<Locator %s>"%str(self.getX())
312
313 def getFunctionSpace(self):
314 """
315 Returns the function space of the Locator.
316 """
317 return self.__function_space
318
319 def getId(self):
320 """
321 Returns the identifier of the location.
322 """
323 return self.__id
324
325 def getX(self):
326 """
327 Returns the exact coordinates of the Locator.
328 """
329 return self(self.getFunctionSpace().getX())
330
331 def __call__(self,data):
332 """
333 Returns the value of data at the Locator of a Data object otherwise
334 the object is returned.
335 """
336 return self.getValue(data)
337
338 def getValue(self,data):
339 """
340 Returns the value of data at the Locator if data is a Data object
341 otherwise the object is returned.
342 """
343 if isinstance(data,escript.Data):
344 if data.getFunctionSpace()==self.getFunctionSpace():
345 #out=data.convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])
346 out=data.convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1],self.getId()[2])
347 else:
348 #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:
351 return out[0]
352 else:
353 return out
354 else:
355 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-7,*args):
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 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,*args):
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 @return: increment dp
424 @rtype: L{escript.Data}
425 @note: this method has to be overwritten by a particular saddle point problem
426 """
427 pass
428
429 def inner(self,p0,p1):
430 """
431 inner product of p0 and p1 approximating p. Typically this returns integrate(p0*p1)
432 @return: inner product of p0 and p1
433 @rtype: C{float}
434 """
435 pass
436
437 def solve(self,u0,p0,tolerance=1.e-6,*args):
438 pass
439 # 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