/[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 1737 - (show annotations)
Fri Aug 29 03:01:29 2008 UTC (10 years, 11 months ago) by gross
File MIME type: text/x-python
File size: 57122 byte(s)
the MaskFromBoundaryTag function didn't make much sense. The argument function_space is replaced by domain.


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 stoppingcriterium2(self,norm_r,norm_b,solver="GMRES",TOL=None):
481 """
482 returns True if the C{norm_r} is C{tolerance}*C{norm_b}
483
484
485 @param norm_r: current residual norm
486 @type norm_r: non-negative C{float}
487 @param norm_b: norm of right hand side
488 @type norm_b: non-negative C{float}
489 @return: C{True} is the stopping criterium is fullfilled. Otherwise C{False} is returned.
490 @rtype: C{bool}
491
492 """
493 if TOL==None:
494 TOL=self.tolerance
495 self.history.append(norm_r)
496 if self.verbose: print "iter: %s: norm(r) = %e"#(len(self.history)-1, self.history[-1])
497 return self.history[-1]<=TOL * norm_b
498
499 def PCG(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
500 """
501 Solver for
502
503 M{Ax=b}
504
505 with a symmetric and positive definite operator A (more details required!).
506 It uses the conjugate gradient method with preconditioner M providing an approximation of A.
507
508 The iteration is terminated if the C{stoppingcriterium} function return C{True}.
509
510 For details on the preconditioned conjugate gradient method see the book:
511
512 Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
513 T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
514 C. Romine, and H. van der Vorst.
515
516 @param b: the right hand side of the liner system. C{b} is altered.
517 @type b: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
518 @param Aprod: returns the value Ax
519 @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}.
520 @param Msolve: solves Mx=r
521 @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
522 type like argument C{x}.
523 @param bilinearform: inner product C{<x,r>}
524 @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}.
525 @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.
526 @type stoppingcriterium: function that returns C{True} or C{False}
527 @param x: an initial guess for the solution. If no C{x} is given 0*b is used.
528 @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
529 @param iter_max: maximum number of iteration steps.
530 @type iter_max: C{int}
531 @return: the solution approximation and the corresponding residual
532 @rtype: C{tuple}
533 @warning: C{b} and C{x} are altered.
534 """
535 iter=0
536 if x==None:
537 x=0*b
538 else:
539 b += (-1)*Aprod(x)
540 r=b
541 rhat=Msolve(r)
542 d = rhat
543 rhat_dot_r = bilinearform(rhat, r)
544 if rhat_dot_r<0: raise NegativeNorm,"negative norm."
545
546 while not stoppingcriterium(math.sqrt(rhat_dot_r),r,x):
547 iter+=1
548 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
549
550 q=Aprod(d)
551 alpha = rhat_dot_r / bilinearform(d, q)
552 x += alpha * d
553 r += (-alpha) * q
554
555 rhat=Msolve(r)
556 rhat_dot_r_new = bilinearform(rhat, r)
557 beta = rhat_dot_r_new / rhat_dot_r
558 rhat+=beta * d
559 d=rhat
560
561 rhat_dot_r = rhat_dot_r_new
562 if rhat_dot_r<0: raise NegativeNorm,"negative norm."
563
564 return x,r
565
566
567 ############################
568 # Added by Artak
569 #################################3
570
571 #Apply a sequence of k Givens rotations, used within gmres codes
572 # vrot=givapp(c, s, vin, k)
573 def givapp(c,s,vin):
574 vrot=vin # warning: vin is altered!!!!
575 if isinstance(c,float):
576 vrot=[c*vrot[0]-s*vrot[1],s*vrot[0]+c*vrot[1]]
577 else:
578 for i in range(len(c)):
579 w1=c[i]*vrot[i]-s[i]*vrot[i+1]
580 w2=s[i]*vrot[i]+c[i]*vrot[i+1]
581 vrot[i:i+2]=w1,w2
582 return vrot
583
584 ##############################################
585 def GMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):
586 ################################################
587 m=iter_restart
588 iter=0
589 xc=x
590 while True:
591 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached"%iter_max
592 xc,stopped=GMRESm(b*1, Aprod, Msolve, bilinearform, stoppingcriterium, x=xc*1, iter_max=iter_max-iter, iter_restart=m)
593 if stopped: break
594 iter+=iter_restart
595 return xc
596
597 def GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):
598 iter=0
599 r=Msolve(b)
600 r_dot_r = bilinearform(r, r)
601 if r_dot_r<0: raise NegativeNorm,"negative norm."
602 norm_b=math.sqrt(r_dot_r)
603
604 if x==None:
605 x=0*b
606 else:
607 r=Msolve(b-Aprod(x))
608 r_dot_r = bilinearform(r, r)
609 if r_dot_r<0: raise NegativeNorm,"negative norm."
610
611 h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)
612 c=numarray.zeros(iter_restart,numarray.Float64)
613 s=numarray.zeros(iter_restart,numarray.Float64)
614 g=numarray.zeros(iter_restart,numarray.Float64)
615 v=[]
616
617 rho=math.sqrt(r_dot_r)
618
619 v.append(r/rho)
620 g[0]=rho
621
622 while not (stoppingcriterium(rho,norm_b) or iter==iter_restart-1):
623
624 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
625
626 p=Msolve(Aprod(v[iter]))
627
628 v.append(p)
629
630 v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
631
632 # Modified Gram-Schmidt
633 for j in range(iter+1):
634 h[j][iter]=bilinearform(v[j],v[iter+1])
635 v[iter+1]-=h[j][iter]*v[j]
636
637 h[iter+1][iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
638 v_norm2=h[iter+1][iter]
639
640 # Reorthogonalize if needed
641 if v_norm1 + 0.001*v_norm2 == v_norm1: #Brown/Hindmarsh condition (default)
642 for j in range(iter+1):
643 hr=bilinearform(v[j],v[iter+1])
644 h[j][iter]=h[j][iter]+hr
645 v[iter+1] -= hr*v[j]
646
647 v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
648 h[iter+1][iter]=v_norm2
649
650 # watch out for happy breakdown
651 if not v_norm2 == 0:
652 v[iter+1]=v[iter+1]/h[iter+1][iter]
653
654 # Form and store the information for the new Givens rotation
655 if iter > 0 :
656 hhat=numarray.zeros(iter+1,numarray.Float64)
657 for i in range(iter+1) : hhat[i]=h[i][iter]
658 hhat=givapp(c[0:iter],s[0:iter],hhat);
659 for i in range(iter+1) : h[i][iter]=hhat[i]
660
661 mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])
662
663 if mu!=0 :
664 c[iter]=h[iter][iter]/mu
665 s[iter]=-h[iter+1][iter]/mu
666 h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter]
667 h[iter+1][iter]=0.0
668 g[iter:iter+2]=givapp(c[iter],s[iter],g[iter:iter+2])
669
670 # Update the residual norm
671
672 rho=abs(g[iter+1])
673 iter+=1
674
675 # At this point either iter > iter_max or rho < tol.
676 # It's time to compute x and leave.
677
678 if iter > 0 :
679 y=numarray.zeros(iter,numarray.Float64)
680 y[iter-1] = g[iter-1] / h[iter-1][iter-1]
681 if iter > 1 :
682 i=iter-2
683 while i>=0 :
684 y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]
685 i=i-1
686 xhat=v[iter-1]*y[iter-1]
687 for i in range(iter-1):
688 xhat += v[i]*y[i]
689 else : xhat=v[0]
690
691 x += xhat
692 if iter<iter_restart-1:
693 stopped=True
694 else:
695 stopped=False
696
697 return x,stopped
698
699
700 ######################################################
701 def dirder(x, w, bilinearform, Aprod, Msolve, f0, b ):
702 ######################################################
703
704 # DIRDER estimates the directional derivative of a function F.
705
706
707 # Hardwired difference increment.
708 #
709 epsnew = 1.0e-07
710 #
711 # Scale the step.
712 #
713 norm_w=math.sqrt(bilinearform(w,w))
714 if norm_w== 0.0:
715 return x/x
716
717 epsnew = epsnew / norm_w
718
719 if norm_w > 0.0:
720 epsnew = epsnew * math.sqrt(bilinearform(x,x))
721 #
722 # DEL and F1 could share the same space if storage
723 # is more important than clarity.
724 #
725
726 DEL = x + epsnew * w
727 f1 = -Msolve(Aprod(DEL))
728 z = ( f1 - f0 ) / epsnew
729 return z
730
731 ######################################################
732 def FDGMRES(f0, Aprod, Msolve, bilinearform, stoppingcriterium, xc=None, x=None, iter_max=100, iter_restart=20,TOL=None):
733 ######################################################
734 b=-f0
735 b_dot_b = bilinearform(b, b)
736 if b_dot_b<0: raise NegativeNorm,"negative norm."
737 norm_b=math.sqrt(b_dot_b)
738
739 r=b
740
741 if x==None:
742 x=0*f0
743 else:
744 r=-dirder(xc,x,bilinearform,Aprod,Msolve,f0,b)-f0
745
746 r_dot_r = bilinearform(r, r)
747 if r_dot_r<0: raise NegativeNorm,"negative norm."
748
749 h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)
750 c=numarray.zeros(iter_restart,numarray.Float64)
751 s=numarray.zeros(iter_restart,numarray.Float64)
752 g=numarray.zeros(iter_restart,numarray.Float64)
753 v=[]
754
755 rho=math.sqrt(r_dot_r)
756
757 v.append(r/rho)
758 g[0]=rho
759 iter=0
760
761 while not (stoppingcriterium(rho,norm_b,solver="FDGMRES",TOL=TOL) or iter==iter_restart-1):
762
763 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
764
765
766 p=dirder(xc, v[iter], bilinearform,Aprod,Msolve,f0,b)
767
768 v.append(p)
769
770 v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
771
772 # Modified Gram-Schmidt
773 for j in range(iter+1):
774 h[j][iter]=bilinearform(v[j],v[iter+1])
775 v[iter+1]+=(-1.)*h[j][iter]*v[j]
776
777 h[iter+1][iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
778 v_norm2=h[iter+1][iter]
779
780
781 # Reorthogonalize if needed
782 if v_norm1 + 0.001*v_norm2 == v_norm1: #Brown/Hindmarsh condition (default)
783 for j in range(iter+1):
784 hr=bilinearform(v[j],v[iter+1])
785 h[j][iter]=h[j][iter]+hr #vhat
786 v[iter+1] +=(-1.)*hr*v[j]
787
788 v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
789 h[iter+1][iter]=v_norm2
790
791 # watch out for happy breakdown
792 if v_norm2 != 0:
793 v[iter+1]=v[iter+1]/h[iter+1][iter]
794
795 # Form and store the information for the new Givens rotation
796 if iter > 0 :
797 hhat=[]
798 for i in range(iter+1) : hhat.append(h[i][iter])
799 hhat=givapp(c[0:iter],s[0:iter],hhat);
800 for i in range(iter+1) : h[i][iter]=hhat[i]
801
802 mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])
803 if mu!=0 :
804 c[iter]=h[iter][iter]/mu
805 s[iter]=-h[iter+1][iter]/mu
806 h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter]
807 h[iter+1][iter]=0.0
808 g[iter:iter+2]=givapp(c[iter],s[iter],g[iter:iter+2])
809
810 # Update the residual norm
811 rho=abs(g[iter+1])
812 iter+=1
813
814 if iter > 0 :
815 y=numarray.zeros(iter,numarray.Float64)
816 y[iter-1] = g[iter-1] / h[iter-1][iter-1]
817 if iter > 1 :
818 i=iter-2
819 while i>=0 :
820 y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]
821 i=i-1
822 xhat=v[iter-1]*y[iter-1]
823 for i in range(iter-1):
824 xhat += v[i]*y[i]
825 else : xhat=v[0]
826
827 x += xhat
828 if iter<iter_restart-1:
829 stopped=True
830 else:
831 stopped=False
832
833 return x,stopped
834
835 #################################################
836 def MINRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
837 #################################################
838 #
839 # minres solves the system of linear equations Ax = b
840 # where A is a symmetric matrix (possibly indefinite or singular)
841 # and b is a given vector.
842 #
843 # "A" may be a dense or sparse matrix (preferably sparse!)
844 # or the name of a function such that
845 # y = A(x)
846 # returns the product y = Ax for any given vector x.
847 #
848 # "M" defines a positive-definite preconditioner M = C C'.
849 # "M" may be a dense or sparse matrix (preferably sparse!)
850 # or the name of a function such that
851 # solves the system My = x for any given vector x.
852 #
853 #
854
855 #------------------------------------------------------------------
856 # Set up y and v for the first Lanczos vector v1.
857 # y = beta1 P' v1, where P = C**(-1).
858 # v is really P' v1.
859 #------------------------------------------------------------------
860 if x==None:
861 x=0*b
862 else:
863 b += (-1)*Aprod(x)
864
865 r1 = b
866 y = Msolve(b)
867 beta1 = bilinearform(y,b)
868
869 if beta1< 0: raise NegativeNorm,"negative norm."
870
871 # If b = 0 exactly, stop with x = 0.
872 if beta1==0: return x*0.
873
874 if beta1> 0:
875 beta1 = math.sqrt(beta1)
876
877 #------------------------------------------------------------------
878 # Initialize quantities.
879 # ------------------------------------------------------------------
880 iter = 0
881 Anorm = 0
882 ynorm = 0
883 oldb = 0
884 beta = beta1
885 dbar = 0
886 epsln = 0
887 phibar = beta1
888 rhs1 = beta1
889 rhs2 = 0
890 rnorm = phibar
891 tnorm2 = 0
892 ynorm2 = 0
893 cs = -1
894 sn = 0
895 w = b*0.
896 w2 = b*0.
897 r2 = r1
898 eps = 0.0001
899
900 #---------------------------------------------------------------------
901 # Main iteration loop.
902 # --------------------------------------------------------------------
903 while not stoppingcriterium(rnorm,Anorm*ynorm,'MINRES'): # checks ||r|| < (||A|| ||x||) * TOL
904
905 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
906 iter = iter + 1
907
908 #-----------------------------------------------------------------
909 # Obtain quantities for the next Lanczos vector vk+1, k = 1, 2,...
910 # The general iteration is similar to the case k = 1 with v0 = 0:
911 #
912 # p1 = Operator * v1 - beta1 * v0,
913 # alpha1 = v1'p1,
914 # q2 = p2 - alpha1 * v1,
915 # beta2^2 = q2'q2,
916 # v2 = (1/beta2) q2.
917 #
918 # Again, y = betak P vk, where P = C**(-1).
919 #-----------------------------------------------------------------
920 s = 1/beta # Normalize previous vector (in y).
921 v = s*y # v = vk if P = I
922
923 y = Aprod(v)
924
925 if iter >= 2:
926 y = y - (beta/oldb)*r1
927
928 alfa = bilinearform(v,y) # alphak
929 y += (- alfa/beta)*r2
930 r1 = r2
931 r2 = y
932 y = Msolve(r2)
933 oldb = beta # oldb = betak
934 beta = bilinearform(y,r2) # beta = betak+1^2
935 if beta < 0: raise NegativeNorm,"negative norm."
936
937 beta = math.sqrt( beta )
938 tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta
939
940 if iter==1: # Initialize a few things.
941 gmax = abs( alfa ) # alpha1
942 gmin = gmax # alpha1
943
944 # Apply previous rotation Qk-1 to get
945 # [deltak epslnk+1] = [cs sn][dbark 0 ]
946 # [gbar k dbar k+1] [sn -cs][alfak betak+1].
947
948 oldeps = epsln
949 delta = cs * dbar + sn * alfa # delta1 = 0 deltak
950 gbar = sn * dbar - cs * alfa # gbar 1 = alfa1 gbar k
951 epsln = sn * beta # epsln2 = 0 epslnk+1
952 dbar = - cs * beta # dbar 2 = beta2 dbar k+1
953
954 # Compute the next plane rotation Qk
955
956 gamma = math.sqrt(gbar*gbar+beta*beta) # gammak
957 gamma = max(gamma,eps)
958 cs = gbar / gamma # ck
959 sn = beta / gamma # sk
960 phi = cs * phibar # phik
961 phibar = sn * phibar # phibark+1
962
963 # Update x.
964
965 denom = 1/gamma
966 w1 = w2
967 w2 = w
968 w = (v - oldeps*w1 - delta*w2) * denom
969 x += phi*w
970
971 # Go round again.
972
973 gmax = max(gmax,gamma)
974 gmin = min(gmin,gamma)
975 z = rhs1 / gamma
976 ynorm2 = z*z + ynorm2
977 rhs1 = rhs2 - delta*z
978 rhs2 = - epsln*z
979
980 # Estimate various norms and test for convergence.
981
982 Anorm = math.sqrt( tnorm2 )
983 ynorm = math.sqrt( ynorm2 )
984
985 rnorm = phibar
986
987 return x
988
989 ######################################
990 def NewtonGMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium,x=None, iter_max=100,iter_restart=20,atol=1.e-2,rtol=1.e-4):
991 #####################################
992 gamma=.9
993 lmaxit=100
994 etamax=.5
995
996 n = 1 #len(x)
997 iter=0
998
999 # evaluate f at the initial iterate
1000 # compute the stop tolerance
1001 #
1002 r=b
1003 if x==None:
1004 x=0*b
1005 else:
1006 r += (-1)*Aprod(x)
1007
1008 f0=-Msolve(r)
1009 fnrm=math.sqrt(bilinearform(f0,f0))/math.sqrt(n)
1010 fnrmo=1
1011 stop_tol=atol + rtol*fnrm
1012 #
1013 # main iteration loop
1014 #
1015 while not stoppingcriterium(fnrm*1,stop_tol,'NewtonGMRES',TOL=1.):
1016
1017 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
1018 #
1019 # keep track of the ratio (rat = fnrm/frnmo)
1020 # of successive residual norms and
1021 # the iteration counter (iter)
1022 #
1023 #rat=fnrm/fnrmo
1024 fnrmo=fnrm
1025 iter+=1
1026 #
1027 # compute the step using a GMRES(m) routine especially designed for this purpose
1028 #
1029 initer=0
1030 while True:
1031 xc,stopped=FDGMRES(f0*1, Aprod, Msolve, bilinearform, stoppingcriterium, xc=x*1, iter_max=lmaxit-initer, iter_restart=iter_restart, TOL=etamax)
1032 if stopped: break
1033 initer+=iter_restart
1034 x+=xc
1035 f0=-Msolve(Aprod(x))
1036 fnrm=math.sqrt(bilinearform(f0,f0))/math.sqrt(n)
1037 rat=fnrm/fnrmo
1038
1039
1040 # adjust eta
1041 #
1042 if etamax > 0:
1043 etaold=etamax
1044 etanew=gamma*rat*rat
1045 if gamma*etaold*etaold > .1 :
1046 etanew=max(etanew,gamma*etaold*etaold)
1047 etamax=min(etanew,etamax)
1048 etamax=max(etamax,.5*stop_tol/fnrm)
1049 return x
1050
1051 def TFQMR(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
1052
1053 # TFQMR solver for linear systems
1054 #
1055 #
1056 # initialization
1057 #
1058 errtol = math.sqrt(bilinearform(b,b))
1059 norm_b=errtol
1060 kmax = iter_max
1061 error = []
1062
1063 if math.sqrt(bilinearform(x,x)) != 0.0:
1064 r = b - Aprod(x)
1065 else:
1066 r = b
1067
1068 r=Msolve(r)
1069
1070 u1=0
1071 u2=0
1072 y1=0
1073 y2=0
1074
1075 w = r
1076 y1 = r
1077 iter = 0
1078 d = 0
1079
1080 v = Msolve(Aprod(y1))
1081 u1 = v
1082
1083 theta = 0.0;
1084 eta = 0.0;
1085 tau = math.sqrt(bilinearform(r,r))
1086 error = [ error, tau ]
1087 rho = tau * tau
1088 m=1
1089 #
1090 # TFQMR iteration
1091 #
1092 # while ( iter < kmax-1 ):
1093
1094 while not stoppingcriterium(tau*math.sqrt ( m + 1 ),norm_b,'TFQMR'):
1095 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
1096
1097 sigma = bilinearform(r,v)
1098
1099 if ( sigma == 0.0 ):
1100 raise 'TFQMR breakdown, sigma=0'
1101
1102
1103 alpha = rho / sigma
1104
1105 for j in range(2):
1106 #
1107 # Compute y2 and u2 only if you have to
1108 #
1109 if ( j == 1 ):
1110 y2 = y1 - alpha * v
1111 u2 = Msolve(Aprod(y2))
1112
1113 m = 2 * (iter+1) - 2 + (j+1)
1114 if j==0:
1115 w = w - alpha * u1
1116 d = y1 + ( theta * theta * eta / alpha ) * d
1117 if j==1:
1118 w = w - alpha * u2
1119 d = y2 + ( theta * theta * eta / alpha ) * d
1120
1121 theta = math.sqrt(bilinearform(w,w))/ tau
1122 c = 1.0 / math.sqrt ( 1.0 + theta * theta )
1123 tau = tau * theta * c
1124 eta = c * c * alpha
1125 x = x + eta * d
1126 #
1127 # Try to terminate the iteration at each pass through the loop
1128 #
1129 # if ( tau * math.sqrt ( m + 1 ) <= errtol ):
1130 # error = [ error, tau ]
1131 # total_iters = iter
1132 # break
1133
1134
1135 if ( rho == 0.0 ):
1136 raise 'TFQMR breakdown, rho=0'
1137
1138
1139 rhon = bilinearform(r,w)
1140 beta = rhon / rho;
1141 rho = rhon;
1142 y1 = w + beta * y2;
1143 u1 = Msolve(Aprod(y1))
1144 v = u1 + beta * ( u2 + beta * v )
1145 error = [ error, tau ]
1146 total_iters = iter
1147
1148 iter = iter + 1
1149
1150 return x
1151
1152
1153 #############################################
1154
1155 class ArithmeticTuple(object):
1156 """
1157 tuple supporting inplace update x+=y and scaling x=a*y where x,y is an ArithmeticTuple and a is a float.
1158
1159 example of usage:
1160
1161 from esys.escript import Data
1162 from numarray import array
1163 a=Data(...)
1164 b=array([1.,4.])
1165 x=ArithmeticTuple(a,b)
1166 y=5.*x
1167
1168 """
1169 def __init__(self,*args):
1170 """
1171 initialize object with elements args.
1172
1173 @param args: tuple of object that support implace add (x+=y) and scaling (x=a*y)
1174 """
1175 self.__items=list(args)
1176
1177 def __len__(self):
1178 """
1179 number of items
1180
1181 @return: number of items
1182 @rtype: C{int}
1183 """
1184 return len(self.__items)
1185
1186 def __getitem__(self,index):
1187 """
1188 get an item
1189
1190 @param index: item to be returned
1191 @type index: C{int}
1192 @return: item with index C{index}
1193 """
1194 return self.__items.__getitem__(index)
1195
1196 def __mul__(self,other):
1197 """
1198 scaling from the right
1199
1200 @param other: scaling factor
1201 @type other: C{float}
1202 @return: itemwise self*other
1203 @rtype: L{ArithmeticTuple}
1204 """
1205 out=[]
1206 other=1.*other
1207 if isinstance(other,float):
1208 for i in range(len(self)):
1209 out.append(self[i]*other)
1210 else:
1211 for i in range(len(self)):
1212 out.append(self[i]*other[i])
1213 return ArithmeticTuple(*tuple(out))
1214
1215 def __rmul__(self,other):
1216 """
1217 scaling from the left
1218
1219 @param other: scaling factor
1220 @type other: C{float}
1221 @return: itemwise other*self
1222 @rtype: L{ArithmeticTuple}
1223 """
1224 out=[]
1225 other=1.*other
1226 if isinstance(other,float):
1227 for i in range(len(self)):
1228 out.append(other*self[i])
1229 else:
1230 for i in range(len(self)):
1231 out.append(other[i]*self[i])
1232 return ArithmeticTuple(*tuple(out))
1233
1234 #########################
1235 # Added by Artak
1236 #########################
1237 def __div__(self,other):
1238 """
1239 dividing from the right
1240
1241 @param other: scaling factor
1242 @type other: C{float}
1243 @return: itemwise self/other
1244 @rtype: L{ArithmeticTuple}
1245 """
1246 out=[]
1247 other=1.*other
1248 if isinstance(other,float):
1249 for i in range(len(self)):
1250 out.append(self[i]*(1./other))
1251 else:
1252 for i in range(len(self)):
1253 out.append(self[i]*(1./other[i]))
1254 return ArithmeticTuple(*tuple(out))
1255
1256 def __rdiv__(self,other):
1257 """
1258 dividing from the left
1259
1260 @param other: scaling factor
1261 @type other: C{float}
1262 @return: itemwise other/self
1263 @rtype: L{ArithmeticTuple}
1264 """
1265 out=[]
1266 other=1.*other
1267 if isinstance(other,float):
1268 for i in range(len(self)):
1269 out.append(other*(1./self[i]))
1270 else:
1271 for i in range(len(self)):
1272 out.append(other[i]*(1./self[i]))
1273 return ArithmeticTuple(*tuple(out))
1274
1275 ##########################################33
1276
1277 def __iadd__(self,other):
1278 """
1279 in-place add of other to self
1280
1281 @param other: increment
1282 @type other: C{ArithmeticTuple}
1283 """
1284 if len(self) != len(other):
1285 raise ValueError,"tuple length must match."
1286 for i in range(len(self)):
1287 self.__items[i]+=other[i]
1288 return self
1289
1290 def __add__(self,other):
1291 """
1292 add other to self
1293
1294 @param other: increment
1295 @type other: C{ArithmeticTuple}
1296 """
1297 # if not isinstance(other,float):
1298 if len(self) != len(other):
1299 raise ValueError,"tuple length must match."
1300 for i in range(len(self)):
1301 self.__items[i]+=other[i]
1302 # else:
1303 # for i in range(len(self)):
1304 # self.__items[i]+=other
1305
1306 return self
1307
1308 def __sub__(self,other):
1309 """
1310 subtract other from self
1311
1312 @param other: increment
1313 @type other: C{ArithmeticTuple}
1314 """
1315 if len(self) != len(other):
1316 raise ValueError,"tuple length must match."
1317 for i in range(len(self)):
1318 self.__items[i]-=other[i]
1319 return self
1320
1321 def __isub__(self,other):
1322 """
1323 subtract other from self
1324
1325 @param other: increment
1326 @type other: C{ArithmeticTuple}
1327 """
1328 if len(self) != len(other):
1329 raise ValueError,"tuple length must match."
1330 for i in range(len(self)):
1331 self.__items[i]-=other[i]
1332 return self
1333
1334 def __neg__(self):
1335 """
1336 negate
1337
1338 """
1339 for i in range(len(self)):
1340 self.__items[i]=-self.__items[i]
1341 return self
1342
1343
1344 class HomogeneousSaddlePointProblem(object):
1345 """
1346 This provides a framwork for solving homogeneous saddle point problem of the form
1347
1348 Av+B^*p=f
1349 Bv =0
1350
1351 for the unknowns v and p and given operators A and B and given right hand side f.
1352 B^* is the adjoint operator of B is the given inner product.
1353
1354 """
1355 def __init__(self,**kwargs):
1356 self.setTolerance()
1357 self.setToleranceReductionFactor()
1358
1359 def initialize(self):
1360 """
1361 initialize the problem (overwrite)
1362 """
1363 pass
1364 def B(self,v):
1365 """
1366 returns Bv (overwrite)
1367 @rtype: equal to the type of p
1368
1369 @note: boundary conditions on p should be zero!
1370 """
1371 pass
1372
1373 def inner(self,p0,p1):
1374 """
1375 returns inner product of two element p0 and p1 (overwrite)
1376
1377 @type p0: equal to the type of p
1378 @type p1: equal to the type of p
1379 @rtype: C{float}
1380
1381 @rtype: equal to the type of p
1382 """
1383 pass
1384
1385 def solve_A(self,u,p):
1386 """
1387 solves Av=f-Au-B^*p with accuracy self.getReducedTolerance() (overwrite)
1388
1389 @rtype: equal to the type of v
1390 @note: boundary conditions on v should be zero!
1391 """
1392 pass
1393
1394 def solve_prec(self,p):
1395 """
1396 provides a preconditioner for BA^{-1}B^* with accuracy self.getReducedTolerance() (overwrite)
1397
1398 @rtype: equal to the type of p
1399 """
1400 pass
1401
1402 def stoppingcriterium(self,Bv,v,p):
1403 """
1404 returns a True if iteration is terminated. (overwrite)
1405
1406 @rtype: C{bool}
1407 """
1408 pass
1409
1410 def __inner(self,p,r):
1411 return self.inner(p,r[1])
1412
1413 def __inner_p(self,p1,p2):
1414 return self.inner(p1,p2)
1415
1416 def __inner_a(self,a1,a2):
1417 return self.inner_a(a1,a2)
1418
1419 def __inner_a1(self,a1,a2):
1420 return self.inner(a1[1],a2[1])
1421
1422 def __stoppingcriterium(self,norm_r,r,p):
1423 return self.stoppingcriterium(r[1],r[0],p)
1424
1425 def __stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None):
1426 return self.stoppingcriterium2(norm_r,norm_b,solver,TOL)
1427
1428 def setTolerance(self,tolerance=1.e-8):
1429 self.__tol=tolerance
1430 def getTolerance(self):
1431 return self.__tol
1432 def setToleranceReductionFactor(self,reduction=0.01):
1433 self.__reduction=reduction
1434 def getSubProblemTolerance(self):
1435 return self.__reduction*self.getTolerance()
1436
1437 def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='PCG',iter_restart=20):
1438 """
1439 solves the saddle point problem using initial guesses v and p.
1440
1441 @param max_iter: maximum number of iteration steps.
1442 """
1443 self.verbose=verbose
1444 self.show_details=show_details and self.verbose
1445
1446 # assume p is known: then v=A^-1(f-B^*p)
1447 # which leads to BA^-1B^*p = BA^-1f
1448
1449 # Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask)
1450 self.__z=v+self.solve_A(v,p*0)
1451 Bz=self.B(self.__z)
1452 #
1453 # solve BA^-1B^*p = Bz
1454 #
1455 #
1456 #
1457 self.iter=0
1458 if solver=='GMRES':
1459 if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter
1460 p=GMRES(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.,iter_restart=iter_restart)
1461 # solve Au=f-B^*p
1462 # A(u-v)=f-B^*p-Av
1463 # u=v+(u-v)
1464 u=v+self.solve_A(v,p)
1465
1466 if solver=='NewtonGMRES':
1467 if self.verbose: print "enter NewtonGMRES method (iter_max=%s)"%max_iter
1468 p=NewtonGMRES(Bz,self.__Aprod_Newton,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.,atol=0,rtol=self.getTolerance())
1469 # solve Au=f-B^*p
1470 # A(u-v)=f-B^*p-Av
1471 # u=v+(u-v)
1472 u=v+self.solve_A(v,p)
1473
1474
1475 if solver=='TFQMR':
1476 if self.verbose: print "enter TFQMR method (iter_max=%s)"%max_iter
1477 p=TFQMR(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)
1478 # solve Au=f-B^*p
1479 # A(u-v)=f-B^*p-Av
1480 # u=v+(u-v)
1481 u=v+self.solve_A(v,p)
1482
1483 if solver=='MINRES':
1484 if self.verbose: print "enter MINRES method (iter_max=%s)"%max_iter
1485 p=MINRES(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)
1486 # solve Au=f-B^*p
1487 # A(u-v)=f-B^*p-Av
1488 # u=v+(u-v)
1489 u=v+self.solve_A(v,p)
1490
1491 if solver=='GMRESC':
1492 if self.verbose: print "enter GMRES coupled method (iter_max=%s)"%max_iter
1493 p0=self.solve_prec1(Bz)
1494 #alfa=(1/self.vol)*util.integrate(util.interpolate(p,escript.Function(self.domain)))
1495 #p-=alfa
1496 x=GMRES(ArithmeticTuple(self.__z*1.,p0*1),self.__Anext,self.__Mempty,self.__inner_a,self.__stoppingcriterium2,iter_max=max_iter, x=ArithmeticTuple(v*1,p*1),iter_restart=20)
1497 #x=NewtonGMRES(ArithmeticTuple(self.__z*1.,p0*1),self.__Aprod_Newton2,self.__Mempty,self.__inner_a,self.__stoppingcriterium2,iter_max=max_iter, x=ArithmeticTuple(v*1,p*1),atol=0,rtol=self.getTolerance())
1498
1499 # solve Au=f-B^*p
1500 # A(u-v)=f-B^*p-Av
1501 # u=v+(u-v)
1502 p=x[1]
1503 u=v+self.solve_A(v,p)
1504 #p=x[1]
1505 #u=x[0]
1506
1507 if solver=='PCG':
1508 # note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv
1509 #
1510 # with Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask)
1511 # A(v-z)= f -Az - B^*p (v-z=0 on fixed_u_mask)
1512 if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter
1513 p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p)
1514 u=r[0]
1515 # print "DIFF=",util.integrate(p)
1516
1517 # print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)
1518
1519 return u,p
1520
1521 def __Msolve(self,r):
1522 return self.solve_prec(r[1])
1523
1524 def __Msolve2(self,r):
1525 return self.solve_prec(r*1)
1526
1527 def __Mempty(self,r):
1528 return r
1529
1530
1531 def __Aprod(self,p):
1532 # return BA^-1B*p
1533 #solve Av =B^*p as Av =f-Az-B^*(-p)
1534 v=self.solve_A(self.__z,-p)
1535 return ArithmeticTuple(v, self.B(v))
1536
1537 def __Anext(self,x):
1538 # return next v,p
1539 #solve Av =-B^*p as Av =f-Az-B^*p
1540
1541 pc=x[1]
1542 v=self.solve_A(self.__z,-pc)
1543 p=self.solve_prec1(self.B(v))
1544
1545 return ArithmeticTuple(v,p)
1546
1547
1548 def __Aprod2(self,p):
1549 # return BA^-1B*p
1550 #solve Av =B^*p as Av =f-Az-B^*(-p)
1551 v=self.solve_A(self.__z,-p)
1552 return self.B(v)
1553
1554 def __Aprod_Newton(self,p):
1555 # return BA^-1B*p - Bz
1556 #solve Av =-B^*p as Av =f-Az-B^*p
1557 v=self.solve_A(self.__z,-p)
1558 return self.B(v-self.__z)
1559
1560 def __Aprod_Newton2(self,x):
1561 # return BA^-1B*p - Bz
1562 #solve Av =-B^*p as Av =f-Az-B^*p
1563 pc=x[1]
1564 v=self.solve_A(self.__z,-pc)
1565 p=self.solve_prec1(self.B(v-self.__z))
1566 return ArithmeticTuple(v,p)
1567
1568 class SaddlePointProblem(object):
1569 """
1570 This implements a solver for a saddlepoint problem
1571
1572 M{f(u,p)=0}
1573 M{g(u)=0}
1574
1575 for u and p. The problem is solved with an inexact Uszawa scheme for p:
1576
1577 M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})}
1578 M{Q_g (p^{k+1}-p^{k}) = g(u^{k+1})}
1579
1580 where Q_f is an approximation of the Jacobiean A_f of f with respect to u and Q_f is an approximation of
1581 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'
1582 Q_g can be difficult, non-linear conjugate gradient method is applied to solve for p, so Q_g plays
1583 in fact the role of a preconditioner.
1584 """
1585 def __init__(self,verbose=False,*args):
1586 """
1587 initializes the problem
1588
1589 @param verbose: switches on the printing out some information
1590 @type verbose: C{bool}
1591 @note: this method may be overwritten by a particular saddle point problem
1592 """
1593 print "SaddlePointProblem should not be used anymore!"
1594 if not isinstance(verbose,bool):
1595 raise TypeError("verbose needs to be of type bool.")
1596 self.__verbose=verbose
1597 self.relaxation=1.
1598
1599 def trace(self,text):
1600 """
1601 prints text if verbose has been set
1602
1603 @param text: a text message
1604 @type text: C{str}
1605 """
1606 if self.__verbose: print "%s: %s"%(str(self),text)
1607
1608 def solve_f(self,u,p,tol=1.e-8):
1609 """
1610 solves
1611
1612 A_f du = f(u,p)
1613
1614 with tolerance C{tol} and return du. A_f is Jacobiean of f with respect to u.
1615
1616 @param u: current approximation of u
1617 @type u: L{escript.Data}
1618 @param p: current approximation of p
1619 @type p: L{escript.Data}
1620 @param tol: tolerance expected for du
1621 @type tol: C{float}
1622 @return: increment du
1623 @rtype: L{escript.Data}
1624 @note: this method has to be overwritten by a particular saddle point problem
1625 """
1626 pass
1627
1628 def solve_g(self,u,tol=1.e-8):
1629 """
1630 solves
1631
1632 Q_g dp = g(u)
1633
1634 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.
1635
1636 @param u: current approximation of u
1637 @type u: L{escript.Data}
1638 @param tol: tolerance expected for dp
1639 @type tol: C{float}
1640 @return: increment dp
1641 @rtype: L{escript.Data}
1642 @note: this method has to be overwritten by a particular saddle point problem
1643 """
1644 pass
1645
1646 def inner(self,p0,p1):
1647 """
1648 inner product of p0 and p1 approximating p. Typically this returns integrate(p0*p1)
1649 @return: inner product of p0 and p1
1650 @rtype: C{float}
1651 """
1652 pass
1653
1654 subiter_max=3
1655 def solve(self,u0,p0,tolerance=1.e-6,tolerance_u=None,iter_max=100,accepted_reduction=0.995,relaxation=None):
1656 """
1657 runs the solver
1658
1659 @param u0: initial guess for C{u}
1660 @type u0: L{esys.escript.Data}
1661 @param p0: initial guess for C{p}
1662 @type p0: L{esys.escript.Data}
1663 @param tolerance: tolerance for relative error in C{u} and C{p}
1664 @type tolerance: positive C{float}
1665 @param tolerance_u: tolerance for relative error in C{u} if different from C{tolerance}
1666 @type tolerance_u: positive C{float}
1667 @param iter_max: maximum number of iteration steps.
1668 @type iter_max: C{int}
1669 @param accepted_reduction: if the norm g cannot be reduced by C{accepted_reduction} backtracking to adapt the
1670 relaxation factor. If C{accepted_reduction=None} no backtracking is used.
1671 @type accepted_reduction: positive C{float} or C{None}
1672 @param relaxation: initial relaxation factor. If C{relaxation==None}, the last relaxation factor is used.
1673 @type relaxation: C{float} or C{None}
1674 """
1675 tol=1.e-2
1676 if tolerance_u==None: tolerance_u=tolerance
1677 if not relaxation==None: self.relaxation=relaxation
1678 if accepted_reduction ==None:
1679 angle_limit=0.
1680 elif accepted_reduction>=1.:
1681 angle_limit=0.
1682 else:
1683 angle_limit=util.sqrt(1-accepted_reduction**2)
1684 self.iter=0
1685 u=u0
1686 p=p0
1687 #
1688 # initialize things:
1689 #
1690 converged=False
1691 #
1692 # start loop:
1693 #
1694 # initial search direction is g
1695 #
1696 while not converged :
1697 if self.iter>iter_max:
1698 raise ArithmeticError("no convergence after %s steps."%self.iter)
1699 f_new=self.solve_f(u,p,tol)
1700 norm_f_new = util.Lsup(f_new)
1701 u_new=u-f_new
1702 g_new=self.solve_g(u_new,tol)
1703 self.iter+=1
1704 norm_g_new = util.sqrt(self.inner(g_new,g_new))
1705 if norm_f_new==0. and norm_g_new==0.: return u, p
1706 if self.iter>1 and not accepted_reduction==None:
1707 #
1708 # did we manage to reduce the norm of G? I
1709 # if not we start a backtracking procedure
1710 #
1711 # print "new/old norm = ",norm_g_new, norm_g, norm_g_new/norm_g
1712 if norm_g_new > accepted_reduction * norm_g:
1713 sub_iter=0
1714 s=self.relaxation
1715 d=g
1716 g_last=g
1717 self.trace(" start substepping: f = %s, g = %s, relaxation = %s."%(norm_f_new, norm_g_new, s))
1718 while sub_iter < self.subiter_max and norm_g_new > accepted_reduction * norm_g:
1719 dg= g_new-g_last
1720 norm_dg=abs(util.sqrt(self.inner(dg,dg))/self.relaxation)
1721 rad=self.inner(g_new,dg)/self.relaxation
1722 # print " ",sub_iter,": rad, norm_dg:",abs(rad), norm_dg*norm_g_new * angle_limit
1723 # print " ",sub_iter,": rad, norm_dg:",rad, norm_dg, norm_g_new, norm_g
1724 if abs(rad) < norm_dg*norm_g_new * angle_limit:
1725 if sub_iter>0: self.trace(" no further improvements expected from backtracking.")
1726 break
1727 r=self.relaxation
1728 self.relaxation= - rad/norm_dg**2
1729 s+=self.relaxation
1730 #####
1731 # a=g_new+self.relaxation*dg/r
1732 # print "predicted new norm = ",util.sqrt(self.inner(a,a)),util.sqrt(self.inner(g_new,g_new)), self.relaxation
1733 #####
1734 g_last=g_new
1735 p+=self.relaxation*d
1736 f_new=self.solve_f(u,p,tol)
1737 u_new=u-f_new
1738 g_new=self.solve_g(u_new,tol)
1739 self.iter+=1
1740 norm_f_new = util.Lsup(f_new)
1741 norm_g_new = util.sqrt(self.inner(g_new,g_new))
1742 # print " ",sub_iter," new g norm",norm_g_new
1743 self.trace(" %s th sub-step: f = %s, g = %s, relaxation = %s."%(sub_iter, norm_f_new, norm_g_new, s))
1744 #
1745 # can we expect reduction of g?
1746 #
1747 # u_last=u_new
1748 sub_iter+=1
1749 self.relaxation=s
1750 #
1751 # check for convergence:
1752 #
1753 norm_u_new = util.Lsup(u_new)
1754 p_new=p+self.relaxation*g_new
1755 norm_p_new = util.sqrt(self.inner(p_new,p_new))
1756 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))
1757
1758 if self.iter>1:
1759 dg2=g_new-g
1760 df2=f_new-f
1761 norm_dg2=util.sqrt(self.inner(dg2,dg2))
1762 norm_df2=util.Lsup(df2)
1763 # print norm_g_new, norm_g, norm_dg, norm_p, tolerance
1764 tol_eq_g=tolerance*norm_dg2/(norm_g*abs(self.relaxation))*norm_p_new
1765 tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new
1766 if norm_g_new <= tol_eq_g and norm_f_new <= tol_eq_f:
1767 converged=True
1768 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
1769 self.trace("convergence after %s steps."%self.iter)
1770 return u,p
1771 # def solve(self,u0,p0,tolerance=1.e-6,iter_max=10,self.relaxation=1.):
1772 # tol=1.e-2
1773 # iter=0
1774 # converged=False
1775 # u=u0*1.
1776 # p=p0*1.
1777 # while not converged and iter<iter_max:
1778 # du=self.solve_f(u,p,tol)
1779 # u-=du
1780 # norm_du=util.Lsup(du)
1781 # norm_u=util.Lsup(u)
1782 #
1783 # dp=self.relaxation*self.solve_g(u,tol)
1784 # p+=dp
1785 # norm_dp=util.sqrt(self.inner(dp,dp))
1786 # norm_p=util.sqrt(self.inner(p,p))
1787 # 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)
1788 # iter+=1
1789 #
1790 # converged = (norm_du <= tolerance*norm_u) and (norm_dp <= tolerance*norm_p)
1791 # if converged:
1792 # print "convergence after %s steps."%iter
1793 # else:
1794 # raise ArithmeticError("no convergence after %s steps."%iter)
1795 #
1796 # return u,p
1797
1798 def MaskFromBoundaryTag(domain,*tags):
1799 """
1800 create a mask on the Solution(domain) function space which one for samples
1801 that touch the boundary tagged by tags.
1802
1803 usage: m=MaskFromBoundaryTag(domain,"left", "right")
1804
1805 @param domain: a given domain
1806 @type domain: L{escript.Domain}
1807 @param tags: boundray tags
1808 @type tags: C{str}
1809 @return: a mask which marks samples that are touching the boundary tagged by any of the given tags.
1810 @rtype: L{escript.Data} of rank 0
1811 """
1812 pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1813 d=escript.Scalar(0.,escript.FunctionOnBoundary(function_space.getDomain()))
1814 for t in tags: d.setTaggedValue(t,1.)
1815 pde.setValue(y=d)
1816 return util.whereNonZero(pde.getRightHandSide())

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26