/[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 1122 - (show annotations)
Tue May 1 03:21:04 2007 UTC (12 years, 3 months ago) by gross
File MIME type: text/x-python
File size: 21624 byte(s)
As the LinearPDE uses coefficients as reduced even if they are handed
over as full the projector runs into a problem when reduced and full
arguments are used in the same projector. Now the projector resets the
coefficients befor starting the projection.


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