/[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 867 - (hide annotations)
Mon Oct 9 06:50:09 2006 UTC (13 years 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 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 gross 867 - SaddlePointProblem - solver for Saddle point problems using the inexact uszawa scheme
11 gross 637
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 jgs 121 """
19    
20 gross 637 __author__="Lutz Gross, l.gross@uq.edu.au"
21 elspeth 609 __copyright__=""" Copyright (c) 2006 by ACcESS MNRF
22     http://www.access.edu.au
23     Primary Business: Queensland, Australia"""
24 elspeth 614 __license__="""Licensed under the Open Software License version 3.0
25     http://www.opensource.org/licenses/osl-3.0.php"""
26 gross 637 __url__="http://www.iservo.edu.au/esys"
27     __version__="$Revision$"
28     __date__="$Date$"
29 elspeth 609
30 gross 637
31 jgs 149 import escript
32     import linearPDEs
33 jgs 121 import numarray
34 jgs 149 import util
35 jgs 121
36 gross 351 class TimeIntegrationManager:
37     """
38     a simple mechanism to manage time dependend values.
39    
40 gross 720 typical usage is::
41 gross 351
42 gross 720 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 gross 351
50 gross 720 @note: currently only p=1 is supported.
51 gross 351 """
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 gross 396 def getValue(self):
72 gross 409 out=self.__v_mem[0]
73     if len(out)==1:
74     return out[0]
75     else:
76     return out
77    
78 gross 351 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 gross 396
113 gross 351
114 jgs 121 class Projector:
115 jgs 149 """
116     The Projector is a factory which projects a discontiuous function onto a
117     continuous function on the a given domain.
118     """
119 jgs 121 def __init__(self, domain, reduce = True, fast=True):
120     """
121 jgs 149 Create a continuous function space projector for a domain.
122 jgs 121
123 jgs 149 @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 jgs 121 """
127 jgs 149 self.__pde = linearPDEs.LinearPDE(domain)
128 jgs 148 if fast:
129 jgs 149 self.__pde.setSolverMethod(linearPDEs.LinearPDE.LUMPING)
130 jgs 121 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 jgs 149 Projects input_data onto a continuous function
141 jgs 121
142 jgs 149 @param input_data: The input_data to be projected.
143 jgs 121 """
144 gross 525 out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())
145 jgs 121 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 gross 525 class NoPDE:
173     """
174     solves the following problem for u:
175 jgs 121
176 gross 525 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 gross 720 @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 gross 525 """
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 gross 720 @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}
208 gross 525 @param Y: right hand side
209 gross 720 @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}
210 gross 525 @param q: location of constraints
211 gross 720 @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}
212 gross 525 @param r: value of solution at locations of constraints
213 gross 720 @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}
214 gross 525 """
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 gross 720 @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}
242 gross 525 @param Y: right hand side
243 gross 720 @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}
244 gross 525 @param q: location of constraints
245 gross 720 @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}
246 gross 525 @param r: value of solution at locations of constraints
247 gross 720 @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}
248 gross 525 """
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 jgs 147 class Locator:
286     """
287 jgs 149 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 jgs 147 """
294    
295     def __init__(self,where,x=numarray.zeros((3,))):
296 jgs 149 """
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 jgs 147 self.__function_space=where
303 jgs 121 else:
304 jgs 149 self.__function_space=escript.ContinuousFunction(where)
305 gross 525 self.__id=util.length(self.__function_space.getX()-x[:self.__function_space.getDim()]).mindp()
306 jgs 121
307 jgs 147 def __str__(self):
308 jgs 149 """
309     Returns the coordinates of the Locator as a string.
310     """
311 jgs 147 return "<Locator %s>"%str(self.getX())
312 jgs 121
313 jgs 147 def getFunctionSpace(self):
314 jgs 149 """
315     Returns the function space of the Locator.
316     """
317 jgs 147 return self.__function_space
318    
319 jgs 121 def getId(self):
320 jgs 149 """
321     Returns the identifier of the location.
322     """
323 jgs 121 return self.__id
324    
325 jgs 147 def getX(self):
326 jgs 149 """
327     Returns the exact coordinates of the Locator.
328     """
329 jgs 147 return self(self.getFunctionSpace().getX())
330 jgs 121
331 jgs 147 def __call__(self,data):
332 jgs 149 """
333     Returns the value of data at the Locator of a Data object otherwise
334     the object is returned.
335     """
336 jgs 147 return self.getValue(data)
337 jgs 121
338 jgs 147 def getValue(self,data):
339 jgs 149 """
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 jgs 147 if data.getFunctionSpace()==self.getFunctionSpace():
345 bcumming 790 #out=data.convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])
346     out=data.convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1],self.getId()[2])
347 jgs 147 else:
348 bcumming 790 #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 jgs 147 if data.getRank()==0:
351     return out[0]
352     else:
353     return out
354     else:
355     return data
356 jgs 149
357 gross 867 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 jgs 149 # 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