/[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 1137 - (show annotations)
Thu May 10 08:11:31 2007 UTC (12 years, 2 months ago) by gross
File MIME type: text/x-python
File size: 21654 byte(s)
This version passes the tests on windows except for 

   * vtk
   * netCDF

The version needs to be tested on altix and linux
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
134 def __call__(self, input_data):
135 """
136 Projects input_data onto a continuous function
137
138 @param input_data: The input_data to be projected.
139 """
140 out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())
141 self.__pde.setValue(Y = escript.Data(), Y_reduced = escript.Data())
142 if input_data.getRank()==0:
143 self.__pde.setValue(Y = input_data)
144 out=self.__pde.getSolution()
145 elif input_data.getRank()==1:
146 for i0 in range(input_data.getShape()[0]):
147 self.__pde.setValue(Y = input_data[i0])
148 out[i0]=self.__pde.getSolution()
149 elif input_data.getRank()==2:
150 for i0 in range(input_data.getShape()[0]):
151 for i1 in range(input_data.getShape()[1]):
152 self.__pde.setValue(Y = input_data[i0,i1])
153 out[i0,i1]=self.__pde.getSolution()
154 elif input_data.getRank()==3:
155 for i0 in range(input_data.getShape()[0]):
156 for i1 in range(input_data.getShape()[1]):
157 for i2 in range(input_data.getShape()[2]):
158 self.__pde.setValue(Y = input_data[i0,i1,i2])
159 out[i0,i1,i2]=self.__pde.getSolution()
160 else:
161 for i0 in range(input_data.getShape()[0]):
162 for i1 in range(input_data.getShape()[1]):
163 for i2 in range(input_data.getShape()[2]):
164 for i3 in range(input_data.getShape()[3]):
165 self.__pde.setValue(Y = input_data[i0,i1,i2,i3])
166 out[i0,i1,i2,i3]=self.__pde.getSolution()
167 return out
168
169 class NoPDE:
170 """
171 solves the following problem for u:
172
173 M{kronecker[i,j]*D[j]*u[j]=Y[i]}
174
175 with constraint
176
177 M{u[j]=r[j]} where M{q[j]>0}
178
179 where D, Y, r and q are given functions of rank 1.
180
181 In the case of scalars this takes the form
182
183 M{D*u=Y}
184
185 with constraint
186
187 M{u=r} where M{q>0}
188
189 where D, Y, r and q are given scalar functions.
190
191 The constraint is overwriting any other condition.
192
193 @note: This class is similar to the L{linearPDEs.LinearPDE} class with A=B=C=X=0 but has the intention
194 that all input parameter are given in L{Solution} or L{ReducedSolution}. The whole
195 thing is a bit strange and I blame Robert.Woodcock@csiro.au for this.
196 """
197 def __init__(self,domain,D=None,Y=None,q=None,r=None):
198 """
199 initialize the problem
200
201 @param domain: domain of the PDE.
202 @type domain: L{Domain}
203 @param D: coefficient of the solution.
204 @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}
205 @param Y: right hand side
206 @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}
207 @param q: location of constraints
208 @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}
209 @param r: value of solution at locations of constraints
210 @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}
211 """
212 self.__domain=domain
213 self.__D=D
214 self.__Y=Y
215 self.__q=q
216 self.__r=r
217 self.__u=None
218 self.__function_space=escript.Solution(self.__domain)
219 def setReducedOn(self):
220 """
221 sets the L{FunctionSpace} of the solution to L{ReducedSolution}
222 """
223 self.__function_space=escript.ReducedSolution(self.__domain)
224 self.__u=None
225
226 def setReducedOff(self):
227 """
228 sets the L{FunctionSpace} of the solution to L{Solution}
229 """
230 self.__function_space=escript.Solution(self.__domain)
231 self.__u=None
232
233 def setValue(self,D=None,Y=None,q=None,r=None):
234 """
235 assigns values to the parameters.
236
237 @param D: coefficient of the solution.
238 @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}
239 @param Y: right hand side
240 @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}
241 @param q: location of constraints
242 @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}
243 @param r: value of solution at locations of constraints
244 @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}
245 """
246 if not D==None:
247 self.__D=D
248 self.__u=None
249 if not Y==None:
250 self.__Y=Y
251 self.__u=None
252 if not q==None:
253 self.__q=q
254 self.__u=None
255 if not r==None:
256 self.__r=r
257 self.__u=None
258
259 def getSolution(self):
260 """
261 returns the solution
262
263 @return: the solution of the problem
264 @rtype: L{Data} object in the L{FunctionSpace} L{Solution} or L{ReducedSolution}.
265 """
266 if self.__u==None:
267 if self.__D==None:
268 raise ValueError,"coefficient D is undefined"
269 D=escript.Data(self.__D,self.__function_space)
270 if D.getRank()>1:
271 raise ValueError,"coefficient D must have rank 0 or 1"
272 if self.__Y==None:
273 self.__u=escript.Data(0.,D.getShape(),self.__function_space)
274 else:
275 self.__u=util.quotient(self.__Y,D)
276 if not self.__q==None:
277 q=util.wherePositive(escript.Data(self.__q,self.__function_space))
278 self.__u*=(1.-q)
279 if not self.__r==None: self.__u+=q*self.__r
280 return self.__u
281
282 class Locator:
283 """
284 Locator provides access to the values of data objects at a given
285 spatial coordinate x.
286
287 In fact, a Locator object finds the sample in the set of samples of a
288 given function space or domain where which is closest to the given
289 point x.
290 """
291
292 def __init__(self,where,x=numarray.zeros((3,))):
293 """
294 Initializes a Locator to access values in Data objects on the Doamin
295 or FunctionSpace where for the sample point which
296 closest to the given point x.
297
298 @param where: function space
299 @type where: L{escript.FunctionSpace}
300 @param x: coefficient of the solution.
301 @type x: L{numarray.NumArray} or C{list} of L{numarray.NumArray}
302 """
303 if isinstance(where,escript.FunctionSpace):
304 self.__function_space=where
305 else:
306 self.__function_space=escript.ContinuousFunction(where)
307 if isinstance(x, list):
308 self.__id=[]
309 for p in x:
310 self.__id.append(util.length(self.__function_space.getX()-p[:self.__function_space.getDim()]).minGlobalDataPoint())
311 else:
312 self.__id=util.length(self.__function_space.getX()-x[:self.__function_space.getDim()]).minGlobalDataPoint()
313
314 def __str__(self):
315 """
316 Returns the coordinates of the Locator as a string.
317 """
318 x=self.getX()
319 if instance(x,list):
320 out="["
321 first=True
322 for xx in x:
323 if not first:
324 out+=","
325 else:
326 first=False
327 out+=str(xx)
328 out+="]>"
329 else:
330 out=str(x)
331 return out
332
333 def getX(self):
334 """
335 Returns the exact coordinates of the Locator.
336 """
337 return self(self.getFunctionSpace().getX())
338
339 def getFunctionSpace(self):
340 """
341 Returns the function space of the Locator.
342 """
343 return self.__function_space
344
345 def getId(self,item=None):
346 """
347 Returns the identifier of the location.
348 """
349 if item == None:
350 return self.__id
351 else:
352 if isinstance(self.__id,list):
353 return self.__id[item]
354 else:
355 return self.__id
356
357
358 def __call__(self,data):
359 """
360 Returns the value of data at the Locator of a Data object otherwise
361 the object is returned.
362 """
363 return self.getValue(data)
364
365 def getValue(self,data):
366 """
367 Returns the value of data at the Locator if data is a Data object
368 otherwise the object is returned.
369 """
370 if isinstance(data,escript.Data):
371 if data.getFunctionSpace()==self.getFunctionSpace():
372 dat=data
373 else:
374 dat=data.interpolate(self.getFunctionSpace())
375 id=self.getId()
376 r=data.getRank()
377 if isinstance(id,list):
378 out=[]
379 for i in id:
380 o=data.getValueOfGlobalDataPoint(*i)
381 if data.getRank()==0:
382 out.append(o[0])
383 else:
384 out.append(o)
385 return out
386 else:
387 out=data.getValueOfGlobalDataPoint(*id)
388 if data.getRank()==0:
389 return out[0]
390 else:
391 return out
392 else:
393 return data
394
395 class SaddlePointProblem(object):
396 """
397 This implements a solver for a saddlepoint problem
398
399 M{f(u,p)=0}
400 M{g(u)=0}
401
402 for u and p. The problem is solved with an inexact Uszawa scheme for p:
403
404 M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})}
405 M{Q_g (p^{k+1}-p^{k}) = g(u^{k+1})}
406
407 where Q_f is an approximation of the Jacobiean A_f of f with respect to u and Q_f is an approximation of
408 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'
409 Q_g can be difficult, non-linear conjugate gradient method is applied to solve for p, so Q_g plays
410 in fact the role of a preconditioner.
411 """
412 def __init__(self,verbose=False,*args):
413 """
414 initializes the problem
415
416 @param verbose: switches on the printing out some information
417 @type verbose: C{bool}
418 @note: this method may be overwritten by a particular saddle point problem
419 """
420 if not isinstance(verbose,bool):
421 raise TypeError("verbose needs to be of type bool.")
422 self.__verbose=verbose
423 self.relaxation=1.
424
425 def trace(self,text):
426 """
427 prints text if verbose has been set
428
429 @param text: a text message
430 @type text: C{str}
431 """
432 if self.__verbose: print "%s: %s"%(str(self),text)
433
434 def solve_f(self,u,p,tol=1.e-8):
435 """
436 solves
437
438 A_f du = f(u,p)
439
440 with tolerance C{tol} and return du. A_f is Jacobiean of f with respect to u.
441
442 @param u: current approximation of u
443 @type u: L{escript.Data}
444 @param p: current approximation of p
445 @type p: L{escript.Data}
446 @param tol: tolerance expected for du
447 @type tol: C{float}
448 @return: increment du
449 @rtype: L{escript.Data}
450 @note: this method has to be overwritten by a particular saddle point problem
451 """
452 pass
453
454 def solve_g(self,u,tol=1.e-8):
455 """
456 solves
457
458 Q_g dp = g(u)
459
460 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.
461
462 @param u: current approximation of u
463 @type u: L{escript.Data}
464 @param tol: tolerance expected for dp
465 @type tol: C{float}
466 @return: increment dp
467 @rtype: L{escript.Data}
468 @note: this method has to be overwritten by a particular saddle point problem
469 """
470 pass
471
472 def inner(self,p0,p1):
473 """
474 inner product of p0 and p1 approximating p. Typically this returns integrate(p0*p1)
475 @return: inner product of p0 and p1
476 @rtype: C{float}
477 """
478 pass
479
480 subiter_max=3
481 def solve(self,u0,p0,tolerance=1.e-6,tolerance_u=None,iter_max=100,accepted_reduction=0.995,relaxation=None):
482 """
483 runs the solver
484
485 @param u0: initial guess for C{u}
486 @type u0: L{esys.escript.Data}
487 @param p0: initial guess for C{p}
488 @type p0: L{esys.escript.Data}
489 @param tolerance: tolerance for relative error in C{u} and C{p}
490 @type tolerance: positive C{float}
491 @param tolerance_u: tolerance for relative error in C{u} if different from C{tolerance}
492 @type tolerance_u: positive C{float}
493 @param iter_max: maximum number of iteration steps.
494 @type iter_max: C{int}
495 @param accepted_reduction: if the norm g cannot be reduced by C{accepted_reduction} backtracking to adapt the
496 relaxation factor. If C{accepted_reduction=None} no backtracking is used.
497 @type accepted_reduction: positive C{float} or C{None}
498 @param relaxation: initial relaxation factor. If C{relaxation==None}, the last relaxation factor is used.
499 @type relaxation: C{float} or C{None}
500 """
501 tol=1.e-2
502 if tolerance_u==None: tolerance_u=tolerance
503 if not relaxation==None: self.relaxation=relaxation
504 if accepted_reduction ==None:
505 angle_limit=0.
506 elif accepted_reduction>=1.:
507 angle_limit=0.
508 else:
509 angle_limit=util.sqrt(1-accepted_reduction**2)
510 self.iter=0
511 u=u0
512 p=p0
513 #
514 # initialize things:
515 #
516 converged=False
517 #
518 # start loop:
519 #
520 # initial search direction is g
521 #
522 while not converged :
523 if self.iter>iter_max:
524 raise ArithmeticError("no convergence after %s steps."%self.iter)
525 f_new=self.solve_f(u,p,tol)
526 norm_f_new = util.Lsup(f_new)
527 u_new=u-f_new
528 g_new=self.solve_g(u_new,tol)
529 self.iter+=1
530 norm_g_new = util.sqrt(self.inner(g_new,g_new))
531 if norm_f_new==0. and norm_g_new==0.: return u, p
532 if self.iter>1 and not accepted_reduction==None:
533 #
534 # did we manage to reduce the norm of G? I
535 # if not we start a backtracking procedure
536 #
537 # print "new/old norm = ",norm_g_new, norm_g, norm_g_new/norm_g
538 if norm_g_new > accepted_reduction * norm_g:
539 sub_iter=0
540 s=self.relaxation
541 d=g
542 g_last=g
543 self.trace(" start substepping: f = %s, g = %s, relaxation = %s."%(norm_f_new, norm_g_new, s))
544 while sub_iter < self.subiter_max and norm_g_new > accepted_reduction * norm_g:
545 dg= g_new-g_last
546 norm_dg=abs(util.sqrt(self.inner(dg,dg))/self.relaxation)
547 rad=self.inner(g_new,dg)/self.relaxation
548 # print " ",sub_iter,": rad, norm_dg:",abs(rad), norm_dg*norm_g_new * angle_limit
549 # print " ",sub_iter,": rad, norm_dg:",rad, norm_dg, norm_g_new, norm_g
550 if abs(rad) < norm_dg*norm_g_new * angle_limit:
551 if sub_iter>0: self.trace(" no further improvements expected from backtracking.")
552 break
553 r=self.relaxation
554 self.relaxation= - rad/norm_dg**2
555 s+=self.relaxation
556 #####
557 # a=g_new+self.relaxation*dg/r
558 # print "predicted new norm = ",util.sqrt(self.inner(a,a)),util.sqrt(self.inner(g_new,g_new)), self.relaxation
559 #####
560 g_last=g_new
561 p+=self.relaxation*d
562 f_new=self.solve_f(u,p,tol)
563 u_new=u-f_new
564 g_new=self.solve_g(u_new,tol)
565 self.iter+=1
566 norm_f_new = util.Lsup(f_new)
567 norm_g_new = util.sqrt(self.inner(g_new,g_new))
568 # print " ",sub_iter," new g norm",norm_g_new
569 self.trace(" %s th sub-step: f = %s, g = %s, relaxation = %s."%(sub_iter, norm_f_new, norm_g_new, s))
570 #
571 # can we expect reduction of g?
572 #
573 # u_last=u_new
574 sub_iter+=1
575 self.relaxation=s
576 #
577 # check for convergence:
578 #
579 norm_u_new = util.Lsup(u_new)
580 p_new=p+self.relaxation*g_new
581 norm_p_new = util.sqrt(self.inner(p_new,p_new))
582 self.trace("%s th step: f = %s, f/u = %s, g = %s, g/p = %s, relaxation = %s."%(self.iter, norm_f_new ,norm_f_new/norm_u_new, norm_g_new, norm_g_new/norm_p_new, self.relaxation))
583
584 if self.iter>1:
585 dg2=g_new-g
586 df2=f_new-f
587 norm_dg2=util.sqrt(self.inner(dg2,dg2))
588 norm_df2=util.Lsup(df2)
589 # print norm_g_new, norm_g, norm_dg, norm_p, tolerance
590 tol_eq_g=tolerance*norm_dg2/(norm_g*abs(self.relaxation))*norm_p_new
591 tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new
592 if norm_g_new <= tol_eq_g and norm_f_new <= tol_eq_f:
593 converged=True
594 f, norm_f, u, norm_u, g, norm_g, p, norm_p = f_new, norm_f_new, u_new, norm_u_new, g_new, norm_g_new, p_new, norm_p_new
595 self.trace("convergence after %s steps."%self.iter)
596 return u,p
597 # def solve(self,u0,p0,tolerance=1.e-6,iter_max=10,self.relaxation=1.):
598 # tol=1.e-2
599 # iter=0
600 # converged=False
601 # u=u0*1.
602 # p=p0*1.
603 # while not converged and iter<iter_max:
604 # du=self.solve_f(u,p,tol)
605 # u-=du
606 # norm_du=util.Lsup(du)
607 # norm_u=util.Lsup(u)
608 #
609 # dp=self.relaxation*self.solve_g(u,tol)
610 # p+=dp
611 # norm_dp=util.sqrt(self.inner(dp,dp))
612 # norm_p=util.sqrt(self.inner(p,p))
613 # print iter,"-th step rel. errror u,p= (%s,%s) (%s,%s)(%s,%s)"%(norm_du,norm_dp,norm_du/norm_u,norm_dp/norm_p,norm_u,norm_p)
614 # iter+=1
615 #
616 # converged = (norm_du <= tolerance*norm_u) and (norm_dp <= tolerance*norm_p)
617 # if converged:
618 # print "convergence after %s steps."%iter
619 # else:
620 # raise ArithmeticError("no convergence after %s steps."%iter)
621 #
622 # return u,p
623
624 # 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