/[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 1467 - (show annotations)
Wed Apr 2 08:10:37 2008 UTC (11 years, 6 months ago) by gross
File MIME type: text/x-python
File size: 38512 byte(s)
a bit of reengineering.
1 #
2 # $Id$
3 #
4 #######################################################
5 #
6 # Copyright 2003-2007 by ACceSS MNRF
7 # Copyright 2007 by University of Queensland
8 #
9 # http://esscc.uq.edu.au
10 # Primary Business: Queensland, Australia
11 # Licensed under the Open Software License version 3.0
12 # http://www.opensource.org/licenses/osl-3.0.php
13 #
14 #######################################################
15 #
16
17 """
18 Provides some tools related to PDEs.
19
20 Currently includes:
21 - Projector - to project a discontinuous
22 - Locator - to trace values in data objects at a certain location
23 - TimeIntegrationManager - to handel extraplotion in time
24 - SaddlePointProblem - solver for Saddle point problems using the inexact uszawa scheme
25
26 @var __author__: name of author
27 @var __copyright__: copyrights
28 @var __license__: licence agreement
29 @var __url__: url entry point on documentation
30 @var __version__: version
31 @var __date__: date of the version
32 """
33
34 __author__="Lutz Gross, l.gross@uq.edu.au"
35 __copyright__=""" Copyright (c) 2006 by ACcESS MNRF
36 http://www.access.edu.au
37 Primary Business: Queensland, Australia"""
38 __license__="""Licensed under the Open Software License version 3.0
39 http://www.opensource.org/licenses/osl-3.0.php"""
40 __url__="http://www.iservo.edu.au/esys"
41 __version__="$Revision$"
42 __date__="$Date$"
43
44
45 import escript
46 import linearPDEs
47 import numarray
48 import util
49 import math
50
51 ##### Added by Artak
52 # from Numeric import zeros,Int,Float64
53 ###################################
54
55
56 class TimeIntegrationManager:
57 """
58 a simple mechanism to manage time dependend values.
59
60 typical usage is::
61
62 dt=0.1 # time increment
63 tm=TimeIntegrationManager(inital_value,p=1)
64 while t<1.
65 v_guess=tm.extrapolate(dt) # extrapolate to t+dt
66 v=...
67 tm.checkin(dt,v)
68 t+=dt
69
70 @note: currently only p=1 is supported.
71 """
72 def __init__(self,*inital_values,**kwargs):
73 """
74 sets up the value manager where inital_value is the initial value and p is order used for extrapolation
75 """
76 if kwargs.has_key("p"):
77 self.__p=kwargs["p"]
78 else:
79 self.__p=1
80 if kwargs.has_key("time"):
81 self.__t=kwargs["time"]
82 else:
83 self.__t=0.
84 self.__v_mem=[inital_values]
85 self.__order=0
86 self.__dt_mem=[]
87 self.__num_val=len(inital_values)
88
89 def getTime(self):
90 return self.__t
91 def getValue(self):
92 out=self.__v_mem[0]
93 if len(out)==1:
94 return out[0]
95 else:
96 return out
97
98 def checkin(self,dt,*values):
99 """
100 adds new values to the manager. the p+1 last value get lost
101 """
102 o=min(self.__order+1,self.__p)
103 self.__order=min(self.__order+1,self.__p)
104 v_mem_new=[values]
105 dt_mem_new=[dt]
106 for i in range(o-1):
107 v_mem_new.append(self.__v_mem[i])
108 dt_mem_new.append(self.__dt_mem[i])
109 v_mem_new.append(self.__v_mem[o-1])
110 self.__order=o
111 self.__v_mem=v_mem_new
112 self.__dt_mem=dt_mem_new
113 self.__t+=dt
114
115 def extrapolate(self,dt):
116 """
117 extrapolates to dt forward in time.
118 """
119 if self.__order==0:
120 out=self.__v_mem[0]
121 else:
122 out=[]
123 for i in range(self.__num_val):
124 out.append((1.+dt/self.__dt_mem[0])*self.__v_mem[0][i]-dt/self.__dt_mem[0]*self.__v_mem[1][i])
125
126 if len(out)==0:
127 return None
128 elif len(out)==1:
129 return out[0]
130 else:
131 return out
132
133
134 class Projector:
135 """
136 The Projector is a factory which projects a discontiuous function onto a
137 continuous function on the a given domain.
138 """
139 def __init__(self, domain, reduce = True, fast=True):
140 """
141 Create a continuous function space projector for a domain.
142
143 @param domain: Domain of the projection.
144 @param reduce: Flag to reduce projection order (default is True)
145 @param fast: Flag to use a fast method based on matrix lumping (default is true)
146 """
147 self.__pde = linearPDEs.LinearPDE(domain)
148 if fast:
149 self.__pde.setSolverMethod(linearPDEs.LinearPDE.LUMPING)
150 self.__pde.setSymmetryOn()
151 self.__pde.setReducedOrderTo(reduce)
152 self.__pde.setValue(D = 1.)
153 return
154
155 def __call__(self, input_data):
156 """
157 Projects input_data onto a continuous function
158
159 @param input_data: The input_data to be projected.
160 """
161 out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())
162 self.__pde.setValue(Y = escript.Data(), Y_reduced = escript.Data())
163 if input_data.getRank()==0:
164 self.__pde.setValue(Y = input_data)
165 out=self.__pde.getSolution()
166 elif input_data.getRank()==1:
167 for i0 in range(input_data.getShape()[0]):
168 self.__pde.setValue(Y = input_data[i0])
169 out[i0]=self.__pde.getSolution()
170 elif input_data.getRank()==2:
171 for i0 in range(input_data.getShape()[0]):
172 for i1 in range(input_data.getShape()[1]):
173 self.__pde.setValue(Y = input_data[i0,i1])
174 out[i0,i1]=self.__pde.getSolution()
175 elif input_data.getRank()==3:
176 for i0 in range(input_data.getShape()[0]):
177 for i1 in range(input_data.getShape()[1]):
178 for i2 in range(input_data.getShape()[2]):
179 self.__pde.setValue(Y = input_data[i0,i1,i2])
180 out[i0,i1,i2]=self.__pde.getSolution()
181 else:
182 for i0 in range(input_data.getShape()[0]):
183 for i1 in range(input_data.getShape()[1]):
184 for i2 in range(input_data.getShape()[2]):
185 for i3 in range(input_data.getShape()[3]):
186 self.__pde.setValue(Y = input_data[i0,i1,i2,i3])
187 out[i0,i1,i2,i3]=self.__pde.getSolution()
188 return out
189
190 class NoPDE:
191 """
192 solves the following problem for u:
193
194 M{kronecker[i,j]*D[j]*u[j]=Y[i]}
195
196 with constraint
197
198 M{u[j]=r[j]} where M{q[j]>0}
199
200 where D, Y, r and q are given functions of rank 1.
201
202 In the case of scalars this takes the form
203
204 M{D*u=Y}
205
206 with constraint
207
208 M{u=r} where M{q>0}
209
210 where D, Y, r and q are given scalar functions.
211
212 The constraint is overwriting any other condition.
213
214 @note: This class is similar to the L{linearPDEs.LinearPDE} class with A=B=C=X=0 but has the intention
215 that all input parameter are given in L{Solution} or L{ReducedSolution}. The whole
216 thing is a bit strange and I blame Robert.Woodcock@csiro.au for this.
217 """
218 def __init__(self,domain,D=None,Y=None,q=None,r=None):
219 """
220 initialize the problem
221
222 @param domain: domain of the PDE.
223 @type domain: L{Domain}
224 @param D: coefficient of the solution.
225 @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}
226 @param Y: right hand side
227 @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}
228 @param q: location of constraints
229 @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}
230 @param r: value of solution at locations of constraints
231 @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}
232 """
233 self.__domain=domain
234 self.__D=D
235 self.__Y=Y
236 self.__q=q
237 self.__r=r
238 self.__u=None
239 self.__function_space=escript.Solution(self.__domain)
240 def setReducedOn(self):
241 """
242 sets the L{FunctionSpace} of the solution to L{ReducedSolution}
243 """
244 self.__function_space=escript.ReducedSolution(self.__domain)
245 self.__u=None
246
247 def setReducedOff(self):
248 """
249 sets the L{FunctionSpace} of the solution to L{Solution}
250 """
251 self.__function_space=escript.Solution(self.__domain)
252 self.__u=None
253
254 def setValue(self,D=None,Y=None,q=None,r=None):
255 """
256 assigns values to the parameters.
257
258 @param D: coefficient of the solution.
259 @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}
260 @param Y: right hand side
261 @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}
262 @param q: location of constraints
263 @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}
264 @param r: value of solution at locations of constraints
265 @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}
266 """
267 if not D==None:
268 self.__D=D
269 self.__u=None
270 if not Y==None:
271 self.__Y=Y
272 self.__u=None
273 if not q==None:
274 self.__q=q
275 self.__u=None
276 if not r==None:
277 self.__r=r
278 self.__u=None
279
280 def getSolution(self):
281 """
282 returns the solution
283
284 @return: the solution of the problem
285 @rtype: L{Data} object in the L{FunctionSpace} L{Solution} or L{ReducedSolution}.
286 """
287 if self.__u==None:
288 if self.__D==None:
289 raise ValueError,"coefficient D is undefined"
290 D=escript.Data(self.__D,self.__function_space)
291 if D.getRank()>1:
292 raise ValueError,"coefficient D must have rank 0 or 1"
293 if self.__Y==None:
294 self.__u=escript.Data(0.,D.getShape(),self.__function_space)
295 else:
296 self.__u=util.quotient(self.__Y,D)
297 if not self.__q==None:
298 q=util.wherePositive(escript.Data(self.__q,self.__function_space))
299 self.__u*=(1.-q)
300 if not self.__r==None: self.__u+=q*self.__r
301 return self.__u
302
303 class Locator:
304 """
305 Locator provides access to the values of data objects at a given
306 spatial coordinate x.
307
308 In fact, a Locator object finds the sample in the set of samples of a
309 given function space or domain where which is closest to the given
310 point x.
311 """
312
313 def __init__(self,where,x=numarray.zeros((3,))):
314 """
315 Initializes a Locator to access values in Data objects on the Doamin
316 or FunctionSpace where for the sample point which
317 closest to the given point x.
318
319 @param where: function space
320 @type where: L{escript.FunctionSpace}
321 @param x: coefficient of the solution.
322 @type x: L{numarray.NumArray} or C{list} of L{numarray.NumArray}
323 """
324 if isinstance(where,escript.FunctionSpace):
325 self.__function_space=where
326 else:
327 self.__function_space=escript.ContinuousFunction(where)
328 if isinstance(x, list):
329 self.__id=[]
330 for p in x:
331 self.__id.append(util.length(self.__function_space.getX()-p[:self.__function_space.getDim()]).minGlobalDataPoint())
332 else:
333 self.__id=util.length(self.__function_space.getX()-x[:self.__function_space.getDim()]).minGlobalDataPoint()
334
335 def __str__(self):
336 """
337 Returns the coordinates of the Locator as a string.
338 """
339 x=self.getX()
340 if instance(x,list):
341 out="["
342 first=True
343 for xx in x:
344 if not first:
345 out+=","
346 else:
347 first=False
348 out+=str(xx)
349 out+="]>"
350 else:
351 out=str(x)
352 return out
353
354 def getX(self):
355 """
356 Returns the exact coordinates of the Locator.
357 """
358 return self(self.getFunctionSpace().getX())
359
360 def getFunctionSpace(self):
361 """
362 Returns the function space of the Locator.
363 """
364 return self.__function_space
365
366 def getId(self,item=None):
367 """
368 Returns the identifier of the location.
369 """
370 if item == None:
371 return self.__id
372 else:
373 if isinstance(self.__id,list):
374 return self.__id[item]
375 else:
376 return self.__id
377
378
379 def __call__(self,data):
380 """
381 Returns the value of data at the Locator of a Data object otherwise
382 the object is returned.
383 """
384 return self.getValue(data)
385
386 def getValue(self,data):
387 """
388 Returns the value of data at the Locator if data is a Data object
389 otherwise the object is returned.
390 """
391 if isinstance(data,escript.Data):
392 if data.getFunctionSpace()==self.getFunctionSpace():
393 dat=data
394 else:
395 dat=data.interpolate(self.getFunctionSpace())
396 id=self.getId()
397 r=data.getRank()
398 if isinstance(id,list):
399 out=[]
400 for i in id:
401 o=data.getValueOfGlobalDataPoint(*i)
402 if data.getRank()==0:
403 out.append(o[0])
404 else:
405 out.append(o)
406 return out
407 else:
408 out=data.getValueOfGlobalDataPoint(*id)
409 if data.getRank()==0:
410 return out[0]
411 else:
412 return out
413 else:
414 return data
415
416 class SolverSchemeException(Exception):
417 """
418 exceptions thrown by solvers
419 """
420 pass
421
422 class IndefinitePreconditioner(SolverSchemeException):
423 """
424 the preconditioner is not positive definite.
425 """
426 pass
427 class MaxIterReached(SolverSchemeException):
428 """
429 maxium number of iteration steps is reached.
430 """
431 pass
432 class IterationBreakDown(SolverSchemeException):
433 """
434 iteration scheme econouters an incurable breakdown.
435 """
436 pass
437 class NegativeNorm(SolverSchemeException):
438 """
439 a norm calculation returns a negative norm.
440 """
441 pass
442
443 class IterationHistory(object):
444 """
445 The IterationHistory class is used to define a stopping criterium. It keeps track of the
446 residual norms. The stoppingcriterium indicates termination if the residual norm has been reduced by
447 a given tolerance.
448 """
449 def __init__(self,tolerance=math.sqrt(util.EPSILON),verbose=False):
450 """
451 Initialization
452
453 @param tolerance: tolerance
454 @type tolerance: positive C{float}
455 @param verbose: switches on the printing out some information
456 @type verbose: C{bool}
457 """
458 if not tolerance>0.:
459 raise ValueError,"tolerance needs to be positive."
460 self.tolerance=tolerance
461 self.verbose=verbose
462 self.history=[]
463 def stoppingcriterium(self,norm_r,r,x):
464 """
465 returns True if the C{norm_r} is C{tolerance}*C{norm_r[0]} where C{norm_r[0]} is the residual norm at the first call.
466
467
468 @param norm_r: current residual norm
469 @type norm_r: non-negative C{float}
470 @param r: current residual (not used)
471 @param x: current solution approximation (not used)
472 @return: C{True} is the stopping criterium is fullfilled. Otherwise C{False} is returned.
473 @rtype: C{bool}
474
475 """
476 self.history.append(norm_r)
477 if self.verbose: print "iter: %s: inner(rhat,r) = %e"%(len(self.history)-1, self.history[-1])
478 return self.history[-1]<=self.tolerance * self.history[0]
479
480 def PCG(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
481 """
482 Solver for
483
484 M{Ax=b}
485
486 with a symmetric and positive definite operator A (more details required!).
487 It uses the conjugate gradient method with preconditioner M providing an approximation of A.
488
489 The iteration is terminated if the C{stoppingcriterium} function return C{True}.
490
491 For details on the preconditioned conjugate gradient method see the book:
492
493 Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
494 T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
495 C. Romine, and H. van der Vorst.
496
497 @param b: the right hand side of the liner system. C{b} is altered.
498 @type b: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
499 @param Aprod: returns the value Ax
500 @type Aprod: function C{Aprod(x)} where C{x} is of the same object like argument C{x}. The returned object needs to be of the same type like argument C{b}.
501 @param Msolve: solves Mx=r
502 @type Msolve: function C{Msolve(r)} where C{r} is of the same type like argument C{b}. The returned object needs to be of the same
503 type like argument C{x}.
504 @param bilinearform: inner product C{<x,r>}
505 @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same type like argument C{x} and C{r} is . The returned value is a C{float}.
506 @param stoppingcriterium: function which returns True if a stopping criterium is meet. C{stoppingcriterium} has the arguments C{norm_r}, C{r} and C{x} giving the current norm of the residual (=C{sqrt(bilinearform(Msolve(r),r)}), the current residual and the current solution approximation. C{stoppingcriterium} is called in each iteration step.
507 @type stoppingcriterium: function that returns C{True} or C{False}
508 @param x: an initial guess for the solution. If no C{x} is given 0*b is used.
509 @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
510 @param iter_max: maximum number of iteration steps.
511 @type iter_max: C{int}
512 @return: the solution approximation and the corresponding residual
513 @rtype: C{tuple}
514 @warning: C{b} and C{x} are altered.
515 """
516 iter=0
517 if x==None:
518 x=0*b
519 else:
520 b += (-1)*Aprod(x)
521 r=b
522 rhat=Msolve(r)
523 d = rhat
524 rhat_dot_r = bilinearform(rhat, r)
525 if rhat_dot_r<0: raise NegativeNorm,"negative norm."
526
527 while not stoppingcriterium(math.sqrt(rhat_dot_r),r,x):
528 iter+=1
529 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
530
531 q=Aprod(d)
532 alpha = rhat_dot_r / bilinearform(d, q)
533 x += alpha * d
534 r += (-alpha) * q
535
536 rhat=Msolve(r)
537 rhat_dot_r_new = bilinearform(rhat, r)
538 beta = rhat_dot_r_new / rhat_dot_r
539 rhat+=beta * d
540 d=rhat
541
542 rhat_dot_r = rhat_dot_r_new
543 if rhat_dot_r<0: raise NegativeNorm,"negative norm."
544
545 return x,r
546
547
548 ############################
549 # Added by Artak
550 #################################3
551
552 #Apply a sequence of k Givens rotations, used within gmres codes
553 # vrot=givapp(c, s, vin, k)
554 def givapp(c,s,vin):
555 vrot=vin # warning: vin is altered!!!!
556 if isinstance(c,float):
557 vrot=[c*vrot[0]-s*vrot[1],s*vrot[0]+c*vrot[1]]
558 else:
559 for i in range(len(c)):
560 w1=c[i]*vrot[i]-s[i]*vrot[i+1]
561 w2=s[i]*vrot[i]+c[i]*vrot[i+1]
562 vrot[i:i+2]=w1,w2
563 return vrot
564
565 def GMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
566 iter=0
567 r=Msolve(b)
568 r_dot_r = bilinearform(r, r)
569 if r_dot_r<0: raise NegativeNorm,"negative norm."
570 norm_b=math.sqrt(r_dot_r)
571
572 if x==None:
573 x=0*b
574 else:
575 r=Msolve(b-Aprod(x))
576 r_dot_r = bilinearform(r, r)
577 if r_dot_r<0: raise NegativeNorm,"negative norm."
578
579 h=numarray.zeros((iter_max,iter_max),numarray.Float64)
580 c=numarray.zeros(iter_max,numarray.Float64)
581 s=numarray.zeros(iter_max,numarray.Float64)
582 g=numarray.zeros(iter_max,numarray.Float64)
583 v=[]
584
585 rho=math.sqrt(r_dot_r)
586 v.append(r/rho)
587 g[0]=rho
588
589 while not stoppingcriterium(rho,norm_b):
590
591 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
592
593
594 p=Msolve(Aprod(v[iter]))
595
596 v.append(p)
597
598 v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
599
600 # Modified Gram-Schmidt
601 for j in range(iter+1):
602 h[j][iter]=bilinearform(v[j],v[iter+1])
603 v[iter+1]+=(-1.)*h[j][iter]*v[j]
604
605 h[iter+1][iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
606 v_norm2=h[iter+1][iter]
607
608
609 # Reorthogonalize if needed
610 if v_norm1 + 0.001*v_norm2 == v_norm1: #Brown/Hindmarsh condition (default)
611 for j in range(iter+1):
612 hr=bilinearform(v[j],v[iter+1])
613 h[j][iter]=h[j][iter]+hr #vhat
614 v[iter+1] +=(-1.)*hr*v[j]
615
616 v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
617 h[iter+1][iter]=v_norm2
618
619 # watch out for happy breakdown
620 if v_norm2 != 0:
621 v[iter+1]=v[iter+1]/h[iter+1][iter]
622
623 # Form and store the information for the new Givens rotation
624 if iter > 0 :
625 hhat=[]
626 for i in range(iter+1) : hhat.append(h[i][iter])
627 hhat=givapp(c[0:iter],s[0:iter],hhat);
628 for i in range(iter+1) : h[i][iter]=hhat[i]
629
630 mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])
631 if mu!=0 :
632 c[iter]=h[iter][iter]/mu
633 s[iter]=-h[iter+1][iter]/mu
634 h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter]
635 h[iter+1][iter]=0.0
636 g[iter:iter+2]=givapp(c[iter],s[iter],g[iter:iter+2])
637
638 # Update the residual norm
639 rho=abs(g[iter+1])
640 iter+=1
641
642 # At this point either iter > iter_max or rho < tol.
643 # It's time to compute x and leave.
644
645 if iter > 0 :
646 y=numarray.zeros(iter,numarray.Float64)
647 y[iter-1] = g[iter-1] / h[iter-1][iter-1]
648 if iter > 1 :
649 i=iter-2
650 while i>=0 :
651 y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]
652 i=i-1
653 xhat=v[iter-1]*y[iter-1]
654 for i in range(iter-1):
655 xhat += v[i]*y[i]
656 else : xhat=v[0]
657
658 x += xhat
659
660 return x
661
662 #############################################
663
664 class ArithmeticTuple(object):
665 """
666 tuple supporting inplace update x+=y and scaling x=a*y where x,y is an ArithmeticTuple and a is a float.
667
668 example of usage:
669
670 from esys.escript import Data
671 from numarray import array
672 a=Data(...)
673 b=array([1.,4.])
674 x=ArithmeticTuple(a,b)
675 y=5.*x
676
677 """
678 def __init__(self,*args):
679 """
680 initialize object with elements args.
681
682 @param args: tuple of object that support implace add (x+=y) and scaling (x=a*y)
683 """
684 self.__items=list(args)
685
686 def __len__(self):
687 """
688 number of items
689
690 @return: number of items
691 @rtype: C{int}
692 """
693 return len(self.__items)
694
695 def __getitem__(self,index):
696 """
697 get an item
698
699 @param index: item to be returned
700 @type index: C{int}
701 @return: item with index C{index}
702 """
703 return self.__items.__getitem__(index)
704
705 def __mul__(self,other):
706 """
707 scaling from the right
708
709 @param other: scaling factor
710 @type other: C{float}
711 @return: itemwise self*other
712 @rtype: L{ArithmeticTuple}
713 """
714 out=[]
715 for i in range(len(self)):
716 out.append(self[i]*other)
717 return ArithmeticTuple(*tuple(out))
718
719 def __rmul__(self,other):
720 """
721 scaling from the left
722
723 @param other: scaling factor
724 @type other: C{float}
725 @return: itemwise other*self
726 @rtype: L{ArithmeticTuple}
727 """
728 out=[]
729 for i in range(len(self)):
730 out.append(other*self[i])
731 return ArithmeticTuple(*tuple(out))
732
733 #########################
734 # Added by Artak
735 #########################
736 def __div__(self,other):
737 """
738 dividing from the right
739
740 @param other: scaling factor
741 @type other: C{float}
742 @return: itemwise self/other
743 @rtype: L{ArithmeticTuple}
744 """
745 out=[]
746 for i in range(len(self)):
747 out.append(self[i]/other)
748 return ArithmeticTuple(*tuple(out))
749
750 def __rdiv__(self,other):
751 """
752 dividing from the left
753
754 @param other: scaling factor
755 @type other: C{float}
756 @return: itemwise other/self
757 @rtype: L{ArithmeticTuple}
758 """
759 out=[]
760 for i in range(len(self)):
761 out.append(other/self[i])
762 return ArithmeticTuple(*tuple(out))
763
764 ##########################################33
765
766 def __iadd__(self,other):
767 """
768 in-place add of other to self
769
770 @param other: increment
771 @type other: C{ArithmeticTuple}
772 """
773 if len(self) != len(other):
774 raise ValueError,"tuple length must match."
775 for i in range(len(self)):
776 self.__items[i]+=other[i]
777 return self
778
779 class HomogeneousSaddlePointProblem(object):
780 """
781 This provides a framwork for solving homogeneous saddle point problem of the form
782
783 Av+B^*p=f
784 Bv =0
785
786 for the unknowns v and p and given operators A and B and given right hand side f.
787 B^* is the adjoint operator of B is the given inner product.
788
789 """
790 def __init__(self,**kwargs):
791 self.setTolerance()
792 self.setToleranceReductionFactor()
793
794 def initialize(self):
795 """
796 initialize the problem (overwrite)
797 """
798 pass
799 def B(self,v):
800 """
801 returns Bv (overwrite)
802 @rtype: equal to the type of p
803
804 @note: boundary conditions on p should be zero!
805 """
806 pass
807
808 def inner(self,p0,p1):
809 """
810 returns inner product of two element p0 and p1 (overwrite)
811
812 @type p0: equal to the type of p
813 @type p1: equal to the type of p
814 @rtype: C{float}
815
816 @rtype: equal to the type of p
817 """
818 pass
819
820 def solve_A(self,u,p):
821 """
822 solves Av=f-Au-B^*p with accuracy self.getReducedTolerance() (overwrite)
823
824 @rtype: equal to the type of v
825 @note: boundary conditions on v should be zero!
826 """
827 pass
828
829 def solve_prec(self,p):
830 """
831 provides a preconditioner for BA^{-1}B^* with accuracy self.getReducedTolerance() (overwrite)
832
833 @rtype: equal to the type of p
834 """
835 pass
836
837 def stoppingcriterium(self,Bv,v,p):
838 """
839 returns a True if iteration is terminated. (overwrite)
840
841 @rtype: C{bool}
842 """
843 pass
844
845 def __inner(self,p,r):
846 return self.inner(p,r[1])
847
848 def __inner_p(self,p1,p2):
849 return self.inner(p1,p2)
850
851 def __stoppingcriterium(self,norm_r,r,p):
852 return self.stoppingcriterium(r[1],r[0],p)
853
854 def __stoppingcriterium_GMRES(self,norm_r,norm_b):
855 return self.stoppingcriterium_GMRES(norm_r,norm_b)
856
857 def setTolerance(self,tolerance=1.e-8):
858 self.__tol=tolerance
859 def getTolerance(self):
860 return self.__tol
861 def setToleranceReductionFactor(self,reduction=0.01):
862 self.__reduction=reduction
863 def getSubProblemTolerance(self):
864 return self.__reduction*self.getTolerance()
865
866 def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='GMRES'):
867 """
868 solves the saddle point problem using initial guesses v and p.
869
870 @param max_iter: maximum number of iteration steps.
871 """
872 self.verbose=verbose
873 self.show_details=show_details and self.verbose
874
875 # Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask)
876
877
878 self.__z=v+self.solve_A(v,p*0)
879
880 Bz=self.B(self.__z)
881 #
882 # solve BA^-1B^*p = Bz
883 #
884 # note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv
885 #
886 # with Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask)
887 # A(v-z)=Az-B^*p-Az = f -Az - B^*p (v-z=0 on fixed_u_mask)
888 #
889 self.iter=0
890 if solver=='GMRES':
891 if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter
892 p=GMRES(Bz,self.__Aprod_GMRES,self.__Msolve_GMRES,self.__inner_p,self.__stoppingcriterium_GMRES,iter_max=max_iter, x=p*1.)
893 # solve Au=f-B^*p
894 # A(u-v)=f-B^*p-Av
895 # u=v+(u-v)
896 u=v+self.solve_A(v,p)
897
898 else:
899 if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter
900 p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p)
901 u=r[0]
902
903 return u,p
904
905 def __Msolve(self,r):
906 return self.solve_prec(r[1])
907
908 def __Msolve_GMRES(self,r):
909 return self.solve_prec(r)
910
911
912 def __Aprod(self,p):
913 # return BA^-1B*p
914 #solve Av =-B^*p as Av =f-Az-B^*p
915 v=self.solve_A(self.__z,p)
916 return ArithmeticTuple(v, -self.B(v))
917
918 def __Aprod_GMRES(self,p):
919 # return BA^-1B*p
920 #solve Av =-B^*p as Av =f-Az-B^*p
921 v=self.solve_A(self.__z,p)
922 return -self.B(v)
923
924 class SaddlePointProblem(object):
925 """
926 This implements a solver for a saddlepoint problem
927
928 M{f(u,p)=0}
929 M{g(u)=0}
930
931 for u and p. The problem is solved with an inexact Uszawa scheme for p:
932
933 M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})}
934 M{Q_g (p^{k+1}-p^{k}) = g(u^{k+1})}
935
936 where Q_f is an approximation of the Jacobiean A_f of f with respect to u and Q_f is an approximation of
937 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'
938 Q_g can be difficult, non-linear conjugate gradient method is applied to solve for p, so Q_g plays
939 in fact the role of a preconditioner.
940 """
941 def __init__(self,verbose=False,*args):
942 """
943 initializes the problem
944
945 @param verbose: switches on the printing out some information
946 @type verbose: C{bool}
947 @note: this method may be overwritten by a particular saddle point problem
948 """
949 if not isinstance(verbose,bool):
950 raise TypeError("verbose needs to be of type bool.")
951 self.__verbose=verbose
952 self.relaxation=1.
953
954 def trace(self,text):
955 """
956 prints text if verbose has been set
957
958 @param text: a text message
959 @type text: C{str}
960 """
961 if self.__verbose: print "%s: %s"%(str(self),text)
962
963 def solve_f(self,u,p,tol=1.e-8):
964 """
965 solves
966
967 A_f du = f(u,p)
968
969 with tolerance C{tol} and return du. A_f is Jacobiean of f with respect to u.
970
971 @param u: current approximation of u
972 @type u: L{escript.Data}
973 @param p: current approximation of p
974 @type p: L{escript.Data}
975 @param tol: tolerance expected for du
976 @type tol: C{float}
977 @return: increment du
978 @rtype: L{escript.Data}
979 @note: this method has to be overwritten by a particular saddle point problem
980 """
981 pass
982
983 def solve_g(self,u,tol=1.e-8):
984 """
985 solves
986
987 Q_g dp = g(u)
988
989 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.
990
991 @param u: current approximation of u
992 @type u: L{escript.Data}
993 @param tol: tolerance expected for dp
994 @type tol: C{float}
995 @return: increment dp
996 @rtype: L{escript.Data}
997 @note: this method has to be overwritten by a particular saddle point problem
998 """
999 pass
1000
1001 def inner(self,p0,p1):
1002 """
1003 inner product of p0 and p1 approximating p. Typically this returns integrate(p0*p1)
1004 @return: inner product of p0 and p1
1005 @rtype: C{float}
1006 """
1007 pass
1008
1009 subiter_max=3
1010 def solve(self,u0,p0,tolerance=1.e-6,tolerance_u=None,iter_max=100,accepted_reduction=0.995,relaxation=None):
1011 """
1012 runs the solver
1013
1014 @param u0: initial guess for C{u}
1015 @type u0: L{esys.escript.Data}
1016 @param p0: initial guess for C{p}
1017 @type p0: L{esys.escript.Data}
1018 @param tolerance: tolerance for relative error in C{u} and C{p}
1019 @type tolerance: positive C{float}
1020 @param tolerance_u: tolerance for relative error in C{u} if different from C{tolerance}
1021 @type tolerance_u: positive C{float}
1022 @param iter_max: maximum number of iteration steps.
1023 @type iter_max: C{int}
1024 @param accepted_reduction: if the norm g cannot be reduced by C{accepted_reduction} backtracking to adapt the
1025 relaxation factor. If C{accepted_reduction=None} no backtracking is used.
1026 @type accepted_reduction: positive C{float} or C{None}
1027 @param relaxation: initial relaxation factor. If C{relaxation==None}, the last relaxation factor is used.
1028 @type relaxation: C{float} or C{None}
1029 """
1030 tol=1.e-2
1031 if tolerance_u==None: tolerance_u=tolerance
1032 if not relaxation==None: self.relaxation=relaxation
1033 if accepted_reduction ==None:
1034 angle_limit=0.
1035 elif accepted_reduction>=1.:
1036 angle_limit=0.
1037 else:
1038 angle_limit=util.sqrt(1-accepted_reduction**2)
1039 self.iter=0
1040 u=u0
1041 p=p0
1042 #
1043 # initialize things:
1044 #
1045 converged=False
1046 #
1047 # start loop:
1048 #
1049 # initial search direction is g
1050 #
1051 while not converged :
1052 if self.iter>iter_max:
1053 raise ArithmeticError("no convergence after %s steps."%self.iter)
1054 f_new=self.solve_f(u,p,tol)
1055 norm_f_new = util.Lsup(f_new)
1056 u_new=u-f_new
1057 g_new=self.solve_g(u_new,tol)
1058 self.iter+=1
1059 norm_g_new = util.sqrt(self.inner(g_new,g_new))
1060 if norm_f_new==0. and norm_g_new==0.: return u, p
1061 if self.iter>1 and not accepted_reduction==None:
1062 #
1063 # did we manage to reduce the norm of G? I
1064 # if not we start a backtracking procedure
1065 #
1066 # print "new/old norm = ",norm_g_new, norm_g, norm_g_new/norm_g
1067 if norm_g_new > accepted_reduction * norm_g:
1068 sub_iter=0
1069 s=self.relaxation
1070 d=g
1071 g_last=g
1072 self.trace(" start substepping: f = %s, g = %s, relaxation = %s."%(norm_f_new, norm_g_new, s))
1073 while sub_iter < self.subiter_max and norm_g_new > accepted_reduction * norm_g:
1074 dg= g_new-g_last
1075 norm_dg=abs(util.sqrt(self.inner(dg,dg))/self.relaxation)
1076 rad=self.inner(g_new,dg)/self.relaxation
1077 # print " ",sub_iter,": rad, norm_dg:",abs(rad), norm_dg*norm_g_new * angle_limit
1078 # print " ",sub_iter,": rad, norm_dg:",rad, norm_dg, norm_g_new, norm_g
1079 if abs(rad) < norm_dg*norm_g_new * angle_limit:
1080 if sub_iter>0: self.trace(" no further improvements expected from backtracking.")
1081 break
1082 r=self.relaxation
1083 self.relaxation= - rad/norm_dg**2
1084 s+=self.relaxation
1085 #####
1086 # a=g_new+self.relaxation*dg/r
1087 # print "predicted new norm = ",util.sqrt(self.inner(a,a)),util.sqrt(self.inner(g_new,g_new)), self.relaxation
1088 #####
1089 g_last=g_new
1090 p+=self.relaxation*d
1091 f_new=self.solve_f(u,p,tol)
1092 u_new=u-f_new
1093 g_new=self.solve_g(u_new,tol)
1094 self.iter+=1
1095 norm_f_new = util.Lsup(f_new)
1096 norm_g_new = util.sqrt(self.inner(g_new,g_new))
1097 # print " ",sub_iter," new g norm",norm_g_new
1098 self.trace(" %s th sub-step: f = %s, g = %s, relaxation = %s."%(sub_iter, norm_f_new, norm_g_new, s))
1099 #
1100 # can we expect reduction of g?
1101 #
1102 # u_last=u_new
1103 sub_iter+=1
1104 self.relaxation=s
1105 #
1106 # check for convergence:
1107 #
1108 norm_u_new = util.Lsup(u_new)
1109 p_new=p+self.relaxation*g_new
1110 norm_p_new = util.sqrt(self.inner(p_new,p_new))
1111 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))
1112
1113 if self.iter>1:
1114 dg2=g_new-g
1115 df2=f_new-f
1116 norm_dg2=util.sqrt(self.inner(dg2,dg2))
1117 norm_df2=util.Lsup(df2)
1118 # print norm_g_new, norm_g, norm_dg, norm_p, tolerance
1119 tol_eq_g=tolerance*norm_dg2/(norm_g*abs(self.relaxation))*norm_p_new
1120 tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new
1121 if norm_g_new <= tol_eq_g and norm_f_new <= tol_eq_f:
1122 converged=True
1123 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
1124 self.trace("convergence after %s steps."%self.iter)
1125 return u,p
1126 # def solve(self,u0,p0,tolerance=1.e-6,iter_max=10,self.relaxation=1.):
1127 # tol=1.e-2
1128 # iter=0
1129 # converged=False
1130 # u=u0*1.
1131 # p=p0*1.
1132 # while not converged and iter<iter_max:
1133 # du=self.solve_f(u,p,tol)
1134 # u-=du
1135 # norm_du=util.Lsup(du)
1136 # norm_u=util.Lsup(u)
1137 #
1138 # dp=self.relaxation*self.solve_g(u,tol)
1139 # p+=dp
1140 # norm_dp=util.sqrt(self.inner(dp,dp))
1141 # norm_p=util.sqrt(self.inner(p,p))
1142 # 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)
1143 # iter+=1
1144 #
1145 # converged = (norm_du <= tolerance*norm_u) and (norm_dp <= tolerance*norm_p)
1146 # if converged:
1147 # print "convergence after %s steps."%iter
1148 # else:
1149 # raise ArithmeticError("no convergence after %s steps."%iter)
1150 #
1151 # return u,p
1152
1153 def MaskFromBoundaryTag(function_space,*tags):
1154 """
1155 create a mask on the given function space which one for samples
1156 that touch the boundary tagged by tags.
1157
1158 usage: m=MaskFromBoundaryTag(Solution(domain),"left", "right")
1159
1160 @param function_space: a given function space
1161 @type function_space: L{escript.FunctionSpace}
1162 @param tags: boundray tags
1163 @type tags: C{str}
1164 @return: a mask which marks samples used by C{function_space} that are touching the
1165 boundary tagged by any of the given tags.
1166 @rtype: L{escript.Data} of rank 0
1167 """
1168 pde=linearPDEs.LinearPDE(function_space.getDomain(),numEquations=1, numSolutions=1)
1169 d=escript.Scalar(0.,escript.FunctionOnBoundary(function_space.getDomain()))
1170 for t in tags: d.setTaggedValue(t,1.)
1171 pde.setValue(y=d)
1172 out=util.whereNonZero(pde.getRightHandSide())
1173 if out.getFunctionSpace() == function_space:
1174 return out
1175 else:
1176 return util.whereNonZero(util.interpolate(out,function_space))
1177
1178
1179

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26