/[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 1809 - (show annotations)
Thu Sep 25 06:43:44 2008 UTC (11 years ago) by ksteube
File MIME type: text/x-python
File size: 57088 byte(s)
Copyright updated in all python files

1
2 ########################################################
3 #
4 # Copyright (c) 2003-2008 by University of Queensland
5 # Earth Systems Science Computational Center (ESSCC)
6 # http://www.uq.edu.au/esscc
7 #
8 # Primary Business: Queensland, Australia
9 # Licensed under the Open Software License version 3.0
10 # http://www.opensource.org/licenses/osl-3.0.php
11 #
12 ########################################################
13
14 __copyright__="""Copyright (c) 2003-2008 by University of Queensland
15 Earth Systems Science Computational Center (ESSCC)
16 http://www.uq.edu.au/esscc
17 Primary Business: Queensland, Australia"""
18 __license__="""Licensed under the Open Software License version 3.0
19 http://www.opensource.org/licenses/osl-3.0.php"""
20 __url__="http://www.uq.edu.au/esscc/escript-finley"
21
22 """
23 Provides some tools related to PDEs.
24
25 Currently includes:
26 - Projector - to project a discontinuous
27 - Locator - to trace values in data objects at a certain location
28 - TimeIntegrationManager - to handel extraplotion in time
29 - SaddlePointProblem - solver for Saddle point problems using the inexact uszawa scheme
30
31 @var __author__: name of author
32 @var __copyright__: copyrights
33 @var __license__: licence agreement
34 @var __url__: url entry point on documentation
35 @var __version__: version
36 @var __date__: date of the version
37 """
38
39 __author__="Lutz Gross, l.gross@uq.edu.au"
40
41
42 import escript
43 import linearPDEs
44 import numarray
45 import util
46 import math
47
48 ##### Added by Artak
49 # from Numeric import zeros,Int,Float64
50 ###################################
51
52
53 class TimeIntegrationManager:
54 """
55 a simple mechanism to manage time dependend values.
56
57 typical usage is::
58
59 dt=0.1 # time increment
60 tm=TimeIntegrationManager(inital_value,p=1)
61 while t<1.
62 v_guess=tm.extrapolate(dt) # extrapolate to t+dt
63 v=...
64 tm.checkin(dt,v)
65 t+=dt
66
67 @note: currently only p=1 is supported.
68 """
69 def __init__(self,*inital_values,**kwargs):
70 """
71 sets up the value manager where inital_value is the initial value and p is order used for extrapolation
72 """
73 if kwargs.has_key("p"):
74 self.__p=kwargs["p"]
75 else:
76 self.__p=1
77 if kwargs.has_key("time"):
78 self.__t=kwargs["time"]
79 else:
80 self.__t=0.
81 self.__v_mem=[inital_values]
82 self.__order=0
83 self.__dt_mem=[]
84 self.__num_val=len(inital_values)
85
86 def getTime(self):
87 return self.__t
88 def getValue(self):
89 out=self.__v_mem[0]
90 if len(out)==1:
91 return out[0]
92 else:
93 return out
94
95 def checkin(self,dt,*values):
96 """
97 adds new values to the manager. the p+1 last value get lost
98 """
99 o=min(self.__order+1,self.__p)
100 self.__order=min(self.__order+1,self.__p)
101 v_mem_new=[values]
102 dt_mem_new=[dt]
103 for i in range(o-1):
104 v_mem_new.append(self.__v_mem[i])
105 dt_mem_new.append(self.__dt_mem[i])
106 v_mem_new.append(self.__v_mem[o-1])
107 self.__order=o
108 self.__v_mem=v_mem_new
109 self.__dt_mem=dt_mem_new
110 self.__t+=dt
111
112 def extrapolate(self,dt):
113 """
114 extrapolates to dt forward in time.
115 """
116 if self.__order==0:
117 out=self.__v_mem[0]
118 else:
119 out=[]
120 for i in range(self.__num_val):
121 out.append((1.+dt/self.__dt_mem[0])*self.__v_mem[0][i]-dt/self.__dt_mem[0]*self.__v_mem[1][i])
122
123 if len(out)==0:
124 return None
125 elif len(out)==1:
126 return out[0]
127 else:
128 return out
129
130
131 class Projector:
132 """
133 The Projector is a factory which projects a discontiuous function onto a
134 continuous function on the a given domain.
135 """
136 def __init__(self, domain, reduce = True, fast=True):
137 """
138 Create a continuous function space projector for a domain.
139
140 @param domain: Domain of the projection.
141 @param reduce: Flag to reduce projection order (default is True)
142 @param fast: Flag to use a fast method based on matrix lumping (default is true)
143 """
144 self.__pde = linearPDEs.LinearPDE(domain)
145 if fast:
146 self.__pde.setSolverMethod(linearPDEs.LinearPDE.LUMPING)
147 self.__pde.setSymmetryOn()
148 self.__pde.setReducedOrderTo(reduce)
149 self.__pde.setValue(D = 1.)
150 return
151
152 def __call__(self, input_data):
153 """
154 Projects input_data onto a continuous function
155
156 @param input_data: The input_data to be projected.
157 """
158 out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())
159 self.__pde.setValue(Y = escript.Data(), Y_reduced = escript.Data())
160 if input_data.getRank()==0:
161 self.__pde.setValue(Y = input_data)
162 out=self.__pde.getSolution()
163 elif input_data.getRank()==1:
164 for i0 in range(input_data.getShape()[0]):
165 self.__pde.setValue(Y = input_data[i0])
166 out[i0]=self.__pde.getSolution()
167 elif input_data.getRank()==2:
168 for i0 in range(input_data.getShape()[0]):
169 for i1 in range(input_data.getShape()[1]):
170 self.__pde.setValue(Y = input_data[i0,i1])
171 out[i0,i1]=self.__pde.getSolution()
172 elif input_data.getRank()==3:
173 for i0 in range(input_data.getShape()[0]):
174 for i1 in range(input_data.getShape()[1]):
175 for i2 in range(input_data.getShape()[2]):
176 self.__pde.setValue(Y = input_data[i0,i1,i2])
177 out[i0,i1,i2]=self.__pde.getSolution()
178 else:
179 for i0 in range(input_data.getShape()[0]):
180 for i1 in range(input_data.getShape()[1]):
181 for i2 in range(input_data.getShape()[2]):
182 for i3 in range(input_data.getShape()[3]):
183 self.__pde.setValue(Y = input_data[i0,i1,i2,i3])
184 out[i0,i1,i2,i3]=self.__pde.getSolution()
185 return out
186
187 class NoPDE:
188 """
189 solves the following problem for u:
190
191 M{kronecker[i,j]*D[j]*u[j]=Y[i]}
192
193 with constraint
194
195 M{u[j]=r[j]} where M{q[j]>0}
196
197 where D, Y, r and q are given functions of rank 1.
198
199 In the case of scalars this takes the form
200
201 M{D*u=Y}
202
203 with constraint
204
205 M{u=r} where M{q>0}
206
207 where D, Y, r and q are given scalar functions.
208
209 The constraint is overwriting any other condition.
210
211 @note: This class is similar to the L{linearPDEs.LinearPDE} class with A=B=C=X=0 but has the intention
212 that all input parameter are given in L{Solution} or L{ReducedSolution}. The whole
213 thing is a bit strange and I blame Robert.Woodcock@csiro.au for this.
214 """
215 def __init__(self,domain,D=None,Y=None,q=None,r=None):
216 """
217 initialize the problem
218
219 @param domain: domain of the PDE.
220 @type domain: L{Domain}
221 @param D: coefficient of the solution.
222 @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}
223 @param Y: right hand side
224 @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}
225 @param q: location of constraints
226 @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}
227 @param r: value of solution at locations of constraints
228 @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}
229 """
230 self.__domain=domain
231 self.__D=D
232 self.__Y=Y
233 self.__q=q
234 self.__r=r
235 self.__u=None
236 self.__function_space=escript.Solution(self.__domain)
237 def setReducedOn(self):
238 """
239 sets the L{FunctionSpace} of the solution to L{ReducedSolution}
240 """
241 self.__function_space=escript.ReducedSolution(self.__domain)
242 self.__u=None
243
244 def setReducedOff(self):
245 """
246 sets the L{FunctionSpace} of the solution to L{Solution}
247 """
248 self.__function_space=escript.Solution(self.__domain)
249 self.__u=None
250
251 def setValue(self,D=None,Y=None,q=None,r=None):
252 """
253 assigns values to the parameters.
254
255 @param D: coefficient of the solution.
256 @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}
257 @param Y: right hand side
258 @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}
259 @param q: location of constraints
260 @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}
261 @param r: value of solution at locations of constraints
262 @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}
263 """
264 if not D==None:
265 self.__D=D
266 self.__u=None
267 if not Y==None:
268 self.__Y=Y
269 self.__u=None
270 if not q==None:
271 self.__q=q
272 self.__u=None
273 if not r==None:
274 self.__r=r
275 self.__u=None
276
277 def getSolution(self):
278 """
279 returns the solution
280
281 @return: the solution of the problem
282 @rtype: L{Data} object in the L{FunctionSpace} L{Solution} or L{ReducedSolution}.
283 """
284 if self.__u==None:
285 if self.__D==None:
286 raise ValueError,"coefficient D is undefined"
287 D=escript.Data(self.__D,self.__function_space)
288 if D.getRank()>1:
289 raise ValueError,"coefficient D must have rank 0 or 1"
290 if self.__Y==None:
291 self.__u=escript.Data(0.,D.getShape(),self.__function_space)
292 else:
293 self.__u=util.quotient(self.__Y,D)
294 if not self.__q==None:
295 q=util.wherePositive(escript.Data(self.__q,self.__function_space))
296 self.__u*=(1.-q)
297 if not self.__r==None: self.__u+=q*self.__r
298 return self.__u
299
300 class Locator:
301 """
302 Locator provides access to the values of data objects at a given
303 spatial coordinate x.
304
305 In fact, a Locator object finds the sample in the set of samples of a
306 given function space or domain where which is closest to the given
307 point x.
308 """
309
310 def __init__(self,where,x=numarray.zeros((3,))):
311 """
312 Initializes a Locator to access values in Data objects on the Doamin
313 or FunctionSpace where for the sample point which
314 closest to the given point x.
315
316 @param where: function space
317 @type where: L{escript.FunctionSpace}
318 @param x: coefficient of the solution.
319 @type x: L{numarray.NumArray} or C{list} of L{numarray.NumArray}
320 """
321 if isinstance(where,escript.FunctionSpace):
322 self.__function_space=where
323 else:
324 self.__function_space=escript.ContinuousFunction(where)
325 if isinstance(x, list):
326 self.__id=[]
327 for p in x:
328 self.__id.append(util.length(self.__function_space.getX()-p[:self.__function_space.getDim()]).minGlobalDataPoint())
329 else:
330 self.__id=util.length(self.__function_space.getX()-x[:self.__function_space.getDim()]).minGlobalDataPoint()
331
332 def __str__(self):
333 """
334 Returns the coordinates of the Locator as a string.
335 """
336 x=self.getX()
337 if instance(x,list):
338 out="["
339 first=True
340 for xx in x:
341 if not first:
342 out+=","
343 else:
344 first=False
345 out+=str(xx)
346 out+="]>"
347 else:
348 out=str(x)
349 return out
350
351 def getX(self):
352 """
353 Returns the exact coordinates of the Locator.
354 """
355 return self(self.getFunctionSpace().getX())
356
357 def getFunctionSpace(self):
358 """
359 Returns the function space of the Locator.
360 """
361 return self.__function_space
362
363 def getId(self,item=None):
364 """
365 Returns the identifier of the location.
366 """
367 if item == None:
368 return self.__id
369 else:
370 if isinstance(self.__id,list):
371 return self.__id[item]
372 else:
373 return self.__id
374
375
376 def __call__(self,data):
377 """
378 Returns the value of data at the Locator of a Data object otherwise
379 the object is returned.
380 """
381 return self.getValue(data)
382
383 def getValue(self,data):
384 """
385 Returns the value of data at the Locator if data is a Data object
386 otherwise the object is returned.
387 """
388 if isinstance(data,escript.Data):
389 if data.getFunctionSpace()==self.getFunctionSpace():
390 dat=data
391 else:
392 dat=data.interpolate(self.getFunctionSpace())
393 id=self.getId()
394 r=data.getRank()
395 if isinstance(id,list):
396 out=[]
397 for i in id:
398 o=data.getValueOfGlobalDataPoint(*i)
399 if data.getRank()==0:
400 out.append(o[0])
401 else:
402 out.append(o)
403 return out
404 else:
405 out=data.getValueOfGlobalDataPoint(*id)
406 if data.getRank()==0:
407 return out[0]
408 else:
409 return out
410 else:
411 return data
412
413 class SolverSchemeException(Exception):
414 """
415 exceptions thrown by solvers
416 """
417 pass
418
419 class IndefinitePreconditioner(SolverSchemeException):
420 """
421 the preconditioner is not positive definite.
422 """
423 pass
424 class MaxIterReached(SolverSchemeException):
425 """
426 maxium number of iteration steps is reached.
427 """
428 pass
429 class IterationBreakDown(SolverSchemeException):
430 """
431 iteration scheme econouters an incurable breakdown.
432 """
433 pass
434 class NegativeNorm(SolverSchemeException):
435 """
436 a norm calculation returns a negative norm.
437 """
438 pass
439
440 class IterationHistory(object):
441 """
442 The IterationHistory class is used to define a stopping criterium. It keeps track of the
443 residual norms. The stoppingcriterium indicates termination if the residual norm has been reduced by
444 a given tolerance.
445 """
446 def __init__(self,tolerance=math.sqrt(util.EPSILON),verbose=False):
447 """
448 Initialization
449
450 @param tolerance: tolerance
451 @type tolerance: positive C{float}
452 @param verbose: switches on the printing out some information
453 @type verbose: C{bool}
454 """
455 if not tolerance>0.:
456 raise ValueError,"tolerance needs to be positive."
457 self.tolerance=tolerance
458 self.verbose=verbose
459 self.history=[]
460 def stoppingcriterium(self,norm_r,r,x):
461 """
462 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.
463
464
465 @param norm_r: current residual norm
466 @type norm_r: non-negative C{float}
467 @param r: current residual (not used)
468 @param x: current solution approximation (not used)
469 @return: C{True} is the stopping criterium is fullfilled. Otherwise C{False} is returned.
470 @rtype: C{bool}
471
472 """
473 self.history.append(norm_r)
474 if self.verbose: print "iter: %s: inner(rhat,r) = %e"#(len(self.history)-1, self.history[-1])
475 return self.history[-1]<=self.tolerance * self.history[0]
476
477 def stoppingcriterium2(self,norm_r,norm_b,solver="GMRES",TOL=None):
478 """
479 returns True if the C{norm_r} is C{tolerance}*C{norm_b}
480
481
482 @param norm_r: current residual norm
483 @type norm_r: non-negative C{float}
484 @param norm_b: norm of right hand side
485 @type norm_b: non-negative C{float}
486 @return: C{True} is the stopping criterium is fullfilled. Otherwise C{False} is returned.
487 @rtype: C{bool}
488
489 """
490 if TOL==None:
491 TOL=self.tolerance
492 self.history.append(norm_r)
493 if self.verbose: print "iter: %s: norm(r) = %e"#(len(self.history)-1, self.history[-1])
494 return self.history[-1]<=TOL * norm_b
495
496 def PCG(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
497 """
498 Solver for
499
500 M{Ax=b}
501
502 with a symmetric and positive definite operator A (more details required!).
503 It uses the conjugate gradient method with preconditioner M providing an approximation of A.
504
505 The iteration is terminated if the C{stoppingcriterium} function return C{True}.
506
507 For details on the preconditioned conjugate gradient method see the book:
508
509 Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
510 T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
511 C. Romine, and H. van der Vorst.
512
513 @param b: the right hand side of the liner system. C{b} is altered.
514 @type b: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
515 @param Aprod: returns the value Ax
516 @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}.
517 @param Msolve: solves Mx=r
518 @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
519 type like argument C{x}.
520 @param bilinearform: inner product C{<x,r>}
521 @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}.
522 @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.
523 @type stoppingcriterium: function that returns C{True} or C{False}
524 @param x: an initial guess for the solution. If no C{x} is given 0*b is used.
525 @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
526 @param iter_max: maximum number of iteration steps.
527 @type iter_max: C{int}
528 @return: the solution approximation and the corresponding residual
529 @rtype: C{tuple}
530 @warning: C{b} and C{x} are altered.
531 """
532 iter=0
533 if x==None:
534 x=0*b
535 else:
536 b += (-1)*Aprod(x)
537 r=b
538 rhat=Msolve(r)
539 d = rhat
540 rhat_dot_r = bilinearform(rhat, r)
541 if rhat_dot_r<0: raise NegativeNorm,"negative norm."
542
543 while not stoppingcriterium(math.sqrt(rhat_dot_r),r,x):
544 iter+=1
545 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
546
547 q=Aprod(d)
548 alpha = rhat_dot_r / bilinearform(d, q)
549 x += alpha * d
550 r += (-alpha) * q
551
552 rhat=Msolve(r)
553 rhat_dot_r_new = bilinearform(rhat, r)
554 beta = rhat_dot_r_new / rhat_dot_r
555 rhat+=beta * d
556 d=rhat
557
558 rhat_dot_r = rhat_dot_r_new
559 if rhat_dot_r<0: raise NegativeNorm,"negative norm."
560
561 return x,r
562
563
564 ############################
565 # Added by Artak
566 #################################3
567
568 #Apply a sequence of k Givens rotations, used within gmres codes
569 # vrot=givapp(c, s, vin, k)
570 def givapp(c,s,vin):
571 vrot=vin # warning: vin is altered!!!!
572 if isinstance(c,float):
573 vrot=[c*vrot[0]-s*vrot[1],s*vrot[0]+c*vrot[1]]
574 else:
575 for i in range(len(c)):
576 w1=c[i]*vrot[i]-s[i]*vrot[i+1]
577 w2=s[i]*vrot[i]+c[i]*vrot[i+1]
578 vrot[i:i+2]=w1,w2
579 return vrot
580
581 ##############################################
582 def GMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):
583 ################################################
584 m=iter_restart
585 iter=0
586 xc=x
587 while True:
588 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached"%iter_max
589 xc,stopped=GMRESm(b*1, Aprod, Msolve, bilinearform, stoppingcriterium, x=xc*1, iter_max=iter_max-iter, iter_restart=m)
590 if stopped: break
591 iter+=iter_restart
592 return xc
593
594 def GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):
595 iter=0
596 r=Msolve(b)
597 r_dot_r = bilinearform(r, r)
598 if r_dot_r<0: raise NegativeNorm,"negative norm."
599 norm_b=math.sqrt(r_dot_r)
600
601 if x==None:
602 x=0*b
603 else:
604 r=Msolve(b-Aprod(x))
605 r_dot_r = bilinearform(r, r)
606 if r_dot_r<0: raise NegativeNorm,"negative norm."
607
608 h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)
609 c=numarray.zeros(iter_restart,numarray.Float64)
610 s=numarray.zeros(iter_restart,numarray.Float64)
611 g=numarray.zeros(iter_restart,numarray.Float64)
612 v=[]
613
614 rho=math.sqrt(r_dot_r)
615
616 v.append(r/rho)
617 g[0]=rho
618
619 while not (stoppingcriterium(rho,norm_b) or iter==iter_restart-1):
620
621 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
622
623 p=Msolve(Aprod(v[iter]))
624
625 v.append(p)
626
627 v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
628
629 # Modified Gram-Schmidt
630 for j in range(iter+1):
631 h[j][iter]=bilinearform(v[j],v[iter+1])
632 v[iter+1]-=h[j][iter]*v[j]
633
634 h[iter+1][iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
635 v_norm2=h[iter+1][iter]
636
637 # Reorthogonalize if needed
638 if v_norm1 + 0.001*v_norm2 == v_norm1: #Brown/Hindmarsh condition (default)
639 for j in range(iter+1):
640 hr=bilinearform(v[j],v[iter+1])
641 h[j][iter]=h[j][iter]+hr
642 v[iter+1] -= hr*v[j]
643
644 v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
645 h[iter+1][iter]=v_norm2
646
647 # watch out for happy breakdown
648 if not v_norm2 == 0:
649 v[iter+1]=v[iter+1]/h[iter+1][iter]
650
651 # Form and store the information for the new Givens rotation
652 if iter > 0 :
653 hhat=numarray.zeros(iter+1,numarray.Float64)
654 for i in range(iter+1) : hhat[i]=h[i][iter]
655 hhat=givapp(c[0:iter],s[0:iter],hhat);
656 for i in range(iter+1) : h[i][iter]=hhat[i]
657
658 mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])
659
660 if mu!=0 :
661 c[iter]=h[iter][iter]/mu
662 s[iter]=-h[iter+1][iter]/mu
663 h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter]
664 h[iter+1][iter]=0.0
665 g[iter:iter+2]=givapp(c[iter],s[iter],g[iter:iter+2])
666
667 # Update the residual norm
668
669 rho=abs(g[iter+1])
670 iter+=1
671
672 # At this point either iter > iter_max or rho < tol.
673 # It's time to compute x and leave.
674
675 if iter > 0 :
676 y=numarray.zeros(iter,numarray.Float64)
677 y[iter-1] = g[iter-1] / h[iter-1][iter-1]
678 if iter > 1 :
679 i=iter-2
680 while i>=0 :
681 y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]
682 i=i-1
683 xhat=v[iter-1]*y[iter-1]
684 for i in range(iter-1):
685 xhat += v[i]*y[i]
686 else : xhat=v[0]
687
688 x += xhat
689 if iter<iter_restart-1:
690 stopped=True
691 else:
692 stopped=False
693
694 return x,stopped
695
696
697 ######################################################
698 def dirder(x, w, bilinearform, Aprod, Msolve, f0, b ):
699 ######################################################
700
701 # DIRDER estimates the directional derivative of a function F.
702
703
704 # Hardwired difference increment.
705 #
706 epsnew = 1.0e-07
707 #
708 # Scale the step.
709 #
710 norm_w=math.sqrt(bilinearform(w,w))
711 if norm_w== 0.0:
712 return x/x
713
714 epsnew = epsnew / norm_w
715
716 if norm_w > 0.0:
717 epsnew = epsnew * math.sqrt(bilinearform(x,x))
718 #
719 # DEL and F1 could share the same space if storage
720 # is more important than clarity.
721 #
722
723 DEL = x + epsnew * w
724 f1 = -Msolve(Aprod(DEL))
725 z = ( f1 - f0 ) / epsnew
726 return z
727
728 ######################################################
729 def FDGMRES(f0, Aprod, Msolve, bilinearform, stoppingcriterium, xc=None, x=None, iter_max=100, iter_restart=20,TOL=None):
730 ######################################################
731 b=-f0
732 b_dot_b = bilinearform(b, b)
733 if b_dot_b<0: raise NegativeNorm,"negative norm."
734 norm_b=math.sqrt(b_dot_b)
735
736 r=b
737
738 if x==None:
739 x=0*f0
740 else:
741 r=-dirder(xc,x,bilinearform,Aprod,Msolve,f0,b)-f0
742
743 r_dot_r = bilinearform(r, r)
744 if r_dot_r<0: raise NegativeNorm,"negative norm."
745
746 h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)
747 c=numarray.zeros(iter_restart,numarray.Float64)
748 s=numarray.zeros(iter_restart,numarray.Float64)
749 g=numarray.zeros(iter_restart,numarray.Float64)
750 v=[]
751
752 rho=math.sqrt(r_dot_r)
753
754 v.append(r/rho)
755 g[0]=rho
756 iter=0
757
758 while not (stoppingcriterium(rho,norm_b,solver="FDGMRES",TOL=TOL) or iter==iter_restart-1):
759
760 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
761
762
763 p=dirder(xc, v[iter], bilinearform,Aprod,Msolve,f0,b)
764
765 v.append(p)
766
767 v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
768
769 # Modified Gram-Schmidt
770 for j in range(iter+1):
771 h[j][iter]=bilinearform(v[j],v[iter+1])
772 v[iter+1]+=(-1.)*h[j][iter]*v[j]
773
774 h[iter+1][iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
775 v_norm2=h[iter+1][iter]
776
777
778 # Reorthogonalize if needed
779 if v_norm1 + 0.001*v_norm2 == v_norm1: #Brown/Hindmarsh condition (default)
780 for j in range(iter+1):
781 hr=bilinearform(v[j],v[iter+1])
782 h[j][iter]=h[j][iter]+hr #vhat
783 v[iter+1] +=(-1.)*hr*v[j]
784
785 v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
786 h[iter+1][iter]=v_norm2
787
788 # watch out for happy breakdown
789 if v_norm2 != 0:
790 v[iter+1]=v[iter+1]/h[iter+1][iter]
791
792 # Form and store the information for the new Givens rotation
793 if iter > 0 :
794 hhat=[]
795 for i in range(iter+1) : hhat.append(h[i][iter])
796 hhat=givapp(c[0:iter],s[0:iter],hhat);
797 for i in range(iter+1) : h[i][iter]=hhat[i]
798
799 mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])
800 if mu!=0 :
801 c[iter]=h[iter][iter]/mu
802 s[iter]=-h[iter+1][iter]/mu
803 h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter]
804 h[iter+1][iter]=0.0
805 g[iter:iter+2]=givapp(c[iter],s[iter],g[iter:iter+2])
806
807 # Update the residual norm
808 rho=abs(g[iter+1])
809 iter+=1
810
811 if iter > 0 :
812 y=numarray.zeros(iter,numarray.Float64)
813 y[iter-1] = g[iter-1] / h[iter-1][iter-1]
814 if iter > 1 :
815 i=iter-2
816 while i>=0 :
817 y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]
818 i=i-1
819 xhat=v[iter-1]*y[iter-1]
820 for i in range(iter-1):
821 xhat += v[i]*y[i]
822 else : xhat=v[0]
823
824 x += xhat
825 if iter<iter_restart-1:
826 stopped=True
827 else:
828 stopped=False
829
830 return x,stopped
831
832 #################################################
833 def MINRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
834 #################################################
835 #
836 # minres solves the system of linear equations Ax = b
837 # where A is a symmetric matrix (possibly indefinite or singular)
838 # and b is a given vector.
839 #
840 # "A" may be a dense or sparse matrix (preferably sparse!)
841 # or the name of a function such that
842 # y = A(x)
843 # returns the product y = Ax for any given vector x.
844 #
845 # "M" defines a positive-definite preconditioner M = C C'.
846 # "M" may be a dense or sparse matrix (preferably sparse!)
847 # or the name of a function such that
848 # solves the system My = x for any given vector x.
849 #
850 #
851
852 #------------------------------------------------------------------
853 # Set up y and v for the first Lanczos vector v1.
854 # y = beta1 P' v1, where P = C**(-1).
855 # v is really P' v1.
856 #------------------------------------------------------------------
857 if x==None:
858 x=0*b
859 else:
860 b += (-1)*Aprod(x)
861
862 r1 = b
863 y = Msolve(b)
864 beta1 = bilinearform(y,b)
865
866 if beta1< 0: raise NegativeNorm,"negative norm."
867
868 # If b = 0 exactly, stop with x = 0.
869 if beta1==0: return x*0.
870
871 if beta1> 0:
872 beta1 = math.sqrt(beta1)
873
874 #------------------------------------------------------------------
875 # Initialize quantities.
876 # ------------------------------------------------------------------
877 iter = 0
878 Anorm = 0
879 ynorm = 0
880 oldb = 0
881 beta = beta1
882 dbar = 0
883 epsln = 0
884 phibar = beta1
885 rhs1 = beta1
886 rhs2 = 0
887 rnorm = phibar
888 tnorm2 = 0
889 ynorm2 = 0
890 cs = -1
891 sn = 0
892 w = b*0.
893 w2 = b*0.
894 r2 = r1
895 eps = 0.0001
896
897 #---------------------------------------------------------------------
898 # Main iteration loop.
899 # --------------------------------------------------------------------
900 while not stoppingcriterium(rnorm,Anorm*ynorm,'MINRES'): # checks ||r|| < (||A|| ||x||) * TOL
901
902 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
903 iter = iter + 1
904
905 #-----------------------------------------------------------------
906 # Obtain quantities for the next Lanczos vector vk+1, k = 1, 2,...
907 # The general iteration is similar to the case k = 1 with v0 = 0:
908 #
909 # p1 = Operator * v1 - beta1 * v0,
910 # alpha1 = v1'p1,
911 # q2 = p2 - alpha1 * v1,
912 # beta2^2 = q2'q2,
913 # v2 = (1/beta2) q2.
914 #
915 # Again, y = betak P vk, where P = C**(-1).
916 #-----------------------------------------------------------------
917 s = 1/beta # Normalize previous vector (in y).
918 v = s*y # v = vk if P = I
919
920 y = Aprod(v)
921
922 if iter >= 2:
923 y = y - (beta/oldb)*r1
924
925 alfa = bilinearform(v,y) # alphak
926 y += (- alfa/beta)*r2
927 r1 = r2
928 r2 = y
929 y = Msolve(r2)
930 oldb = beta # oldb = betak
931 beta = bilinearform(y,r2) # beta = betak+1^2
932 if beta < 0: raise NegativeNorm,"negative norm."
933
934 beta = math.sqrt( beta )
935 tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta
936
937 if iter==1: # Initialize a few things.
938 gmax = abs( alfa ) # alpha1
939 gmin = gmax # alpha1
940
941 # Apply previous rotation Qk-1 to get
942 # [deltak epslnk+1] = [cs sn][dbark 0 ]
943 # [gbar k dbar k+1] [sn -cs][alfak betak+1].
944
945 oldeps = epsln
946 delta = cs * dbar + sn * alfa # delta1 = 0 deltak
947 gbar = sn * dbar - cs * alfa # gbar 1 = alfa1 gbar k
948 epsln = sn * beta # epsln2 = 0 epslnk+1
949 dbar = - cs * beta # dbar 2 = beta2 dbar k+1
950
951 # Compute the next plane rotation Qk
952
953 gamma = math.sqrt(gbar*gbar+beta*beta) # gammak
954 gamma = max(gamma,eps)
955 cs = gbar / gamma # ck
956 sn = beta / gamma # sk
957 phi = cs * phibar # phik
958 phibar = sn * phibar # phibark+1
959
960 # Update x.
961
962 denom = 1/gamma
963 w1 = w2
964 w2 = w
965 w = (v - oldeps*w1 - delta*w2) * denom
966 x += phi*w
967
968 # Go round again.
969
970 gmax = max(gmax,gamma)
971 gmin = min(gmin,gamma)
972 z = rhs1 / gamma
973 ynorm2 = z*z + ynorm2
974 rhs1 = rhs2 - delta*z
975 rhs2 = - epsln*z
976
977 # Estimate various norms and test for convergence.
978
979 Anorm = math.sqrt( tnorm2 )
980 ynorm = math.sqrt( ynorm2 )
981
982 rnorm = phibar
983
984 return x
985
986 ######################################
987 def NewtonGMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium,x=None, iter_max=100,iter_restart=20,atol=1.e-2,rtol=1.e-4):
988 #####################################
989 gamma=.9
990 lmaxit=100
991 etamax=.5
992
993 n = 1 #len(x)
994 iter=0
995
996 # evaluate f at the initial iterate
997 # compute the stop tolerance
998 #
999 r=b
1000 if x==None:
1001 x=0*b
1002 else:
1003 r += (-1)*Aprod(x)
1004
1005 f0=-Msolve(r)
1006 fnrm=math.sqrt(bilinearform(f0,f0))/math.sqrt(n)
1007 fnrmo=1
1008 stop_tol=atol + rtol*fnrm
1009 #
1010 # main iteration loop
1011 #
1012 while not stoppingcriterium(fnrm*1,stop_tol,'NewtonGMRES',TOL=1.):
1013
1014 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
1015 #
1016 # keep track of the ratio (rat = fnrm/frnmo)
1017 # of successive residual norms and
1018 # the iteration counter (iter)
1019 #
1020 #rat=fnrm/fnrmo
1021 fnrmo=fnrm
1022 iter+=1
1023 #
1024 # compute the step using a GMRES(m) routine especially designed for this purpose
1025 #
1026 initer=0
1027 while True:
1028 xc,stopped=FDGMRES(f0*1, Aprod, Msolve, bilinearform, stoppingcriterium, xc=x*1, iter_max=lmaxit-initer, iter_restart=iter_restart, TOL=etamax)
1029 if stopped: break
1030 initer+=iter_restart
1031 x+=xc
1032 f0=-Msolve(Aprod(x))
1033 fnrm=math.sqrt(bilinearform(f0,f0))/math.sqrt(n)
1034 rat=fnrm/fnrmo
1035
1036
1037 # adjust eta
1038 #
1039 if etamax > 0:
1040 etaold=etamax
1041 etanew=gamma*rat*rat
1042 if gamma*etaold*etaold > .1 :
1043 etanew=max(etanew,gamma*etaold*etaold)
1044 etamax=min(etanew,etamax)
1045 etamax=max(etamax,.5*stop_tol/fnrm)
1046 return x
1047
1048 def TFQMR(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
1049
1050 # TFQMR solver for linear systems
1051 #
1052 #
1053 # initialization
1054 #
1055 errtol = math.sqrt(bilinearform(b,b))
1056 norm_b=errtol
1057 kmax = iter_max
1058 error = []
1059
1060 if math.sqrt(bilinearform(x,x)) != 0.0:
1061 r = b - Aprod(x)
1062 else:
1063 r = b
1064
1065 r=Msolve(r)
1066
1067 u1=0
1068 u2=0
1069 y1=0
1070 y2=0
1071
1072 w = r
1073 y1 = r
1074 iter = 0
1075 d = 0
1076
1077 v = Msolve(Aprod(y1))
1078 u1 = v
1079
1080 theta = 0.0;
1081 eta = 0.0;
1082 tau = math.sqrt(bilinearform(r,r))
1083 error = [ error, tau ]
1084 rho = tau * tau
1085 m=1
1086 #
1087 # TFQMR iteration
1088 #
1089 # while ( iter < kmax-1 ):
1090
1091 while not stoppingcriterium(tau*math.sqrt ( m + 1 ),norm_b,'TFQMR'):
1092 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
1093
1094 sigma = bilinearform(r,v)
1095
1096 if ( sigma == 0.0 ):
1097 raise 'TFQMR breakdown, sigma=0'
1098
1099
1100 alpha = rho / sigma
1101
1102 for j in range(2):
1103 #
1104 # Compute y2 and u2 only if you have to
1105 #
1106 if ( j == 1 ):
1107 y2 = y1 - alpha * v
1108 u2 = Msolve(Aprod(y2))
1109
1110 m = 2 * (iter+1) - 2 + (j+1)
1111 if j==0:
1112 w = w - alpha * u1
1113 d = y1 + ( theta * theta * eta / alpha ) * d
1114 if j==1:
1115 w = w - alpha * u2
1116 d = y2 + ( theta * theta * eta / alpha ) * d
1117
1118 theta = math.sqrt(bilinearform(w,w))/ tau
1119 c = 1.0 / math.sqrt ( 1.0 + theta * theta )
1120 tau = tau * theta * c
1121 eta = c * c * alpha
1122 x = x + eta * d
1123 #
1124 # Try to terminate the iteration at each pass through the loop
1125 #
1126 # if ( tau * math.sqrt ( m + 1 ) <= errtol ):
1127 # error = [ error, tau ]
1128 # total_iters = iter
1129 # break
1130
1131
1132 if ( rho == 0.0 ):
1133 raise 'TFQMR breakdown, rho=0'
1134
1135
1136 rhon = bilinearform(r,w)
1137 beta = rhon / rho;
1138 rho = rhon;
1139 y1 = w + beta * y2;
1140 u1 = Msolve(Aprod(y1))
1141 v = u1 + beta * ( u2 + beta * v )
1142 error = [ error, tau ]
1143 total_iters = iter
1144
1145 iter = iter + 1
1146
1147 return x
1148
1149
1150 #############################################
1151
1152 class ArithmeticTuple(object):
1153 """
1154 tuple supporting inplace update x+=y and scaling x=a*y where x,y is an ArithmeticTuple and a is a float.
1155
1156 example of usage:
1157
1158 from esys.escript import Data
1159 from numarray import array
1160 a=Data(...)
1161 b=array([1.,4.])
1162 x=ArithmeticTuple(a,b)
1163 y=5.*x
1164
1165 """
1166 def __init__(self,*args):
1167 """
1168 initialize object with elements args.
1169
1170 @param args: tuple of object that support implace add (x+=y) and scaling (x=a*y)
1171 """
1172 self.__items=list(args)
1173
1174 def __len__(self):
1175 """
1176 number of items
1177
1178 @return: number of items
1179 @rtype: C{int}
1180 """
1181 return len(self.__items)
1182
1183 def __getitem__(self,index):
1184 """
1185 get an item
1186
1187 @param index: item to be returned
1188 @type index: C{int}
1189 @return: item with index C{index}
1190 """
1191 return self.__items.__getitem__(index)
1192
1193 def __mul__(self,other):
1194 """
1195 scaling from the right
1196
1197 @param other: scaling factor
1198 @type other: C{float}
1199 @return: itemwise self*other
1200 @rtype: L{ArithmeticTuple}
1201 """
1202 out=[]
1203 other=1.*other
1204 if isinstance(other,float):
1205 for i in range(len(self)):
1206 out.append(self[i]*other)
1207 else:
1208 for i in range(len(self)):
1209 out.append(self[i]*other[i])
1210 return ArithmeticTuple(*tuple(out))
1211
1212 def __rmul__(self,other):
1213 """
1214 scaling from the left
1215
1216 @param other: scaling factor
1217 @type other: C{float}
1218 @return: itemwise other*self
1219 @rtype: L{ArithmeticTuple}
1220 """
1221 out=[]
1222 other=1.*other
1223 if isinstance(other,float):
1224 for i in range(len(self)):
1225 out.append(other*self[i])
1226 else:
1227 for i in range(len(self)):
1228 out.append(other[i]*self[i])
1229 return ArithmeticTuple(*tuple(out))
1230
1231 #########################
1232 # Added by Artak
1233 #########################
1234 def __div__(self,other):
1235 """
1236 dividing from the right
1237
1238 @param other: scaling factor
1239 @type other: C{float}
1240 @return: itemwise self/other
1241 @rtype: L{ArithmeticTuple}
1242 """
1243 out=[]
1244 other=1.*other
1245 if isinstance(other,float):
1246 for i in range(len(self)):
1247 out.append(self[i]*(1./other))
1248 else:
1249 for i in range(len(self)):
1250 out.append(self[i]*(1./other[i]))
1251 return ArithmeticTuple(*tuple(out))
1252
1253 def __rdiv__(self,other):
1254 """
1255 dividing from the left
1256
1257 @param other: scaling factor
1258 @type other: C{float}
1259 @return: itemwise other/self
1260 @rtype: L{ArithmeticTuple}
1261 """
1262 out=[]
1263 other=1.*other
1264 if isinstance(other,float):
1265 for i in range(len(self)):
1266 out.append(other*(1./self[i]))
1267 else:
1268 for i in range(len(self)):
1269 out.append(other[i]*(1./self[i]))
1270 return ArithmeticTuple(*tuple(out))
1271
1272 ##########################################33
1273
1274 def __iadd__(self,other):
1275 """
1276 in-place add of other to self
1277
1278 @param other: increment
1279 @type other: C{ArithmeticTuple}
1280 """
1281 if len(self) != len(other):
1282 raise ValueError,"tuple length must match."
1283 for i in range(len(self)):
1284 self.__items[i]+=other[i]
1285 return self
1286
1287 def __add__(self,other):
1288 """
1289 add other to self
1290
1291 @param other: increment
1292 @type other: C{ArithmeticTuple}
1293 """
1294 # if not isinstance(other,float):
1295 if len(self) != len(other):
1296 raise ValueError,"tuple length must match."
1297 for i in range(len(self)):
1298 self.__items[i]+=other[i]
1299 # else:
1300 # for i in range(len(self)):
1301 # self.__items[i]+=other
1302
1303 return self
1304
1305 def __sub__(self,other):
1306 """
1307 subtract other from self
1308
1309 @param other: increment
1310 @type other: C{ArithmeticTuple}
1311 """
1312 if len(self) != len(other):
1313 raise ValueError,"tuple length must match."
1314 for i in range(len(self)):
1315 self.__items[i]-=other[i]
1316 return self
1317
1318 def __isub__(self,other):
1319 """
1320 subtract other from self
1321
1322 @param other: increment
1323 @type other: C{ArithmeticTuple}
1324 """
1325 if len(self) != len(other):
1326 raise ValueError,"tuple length must match."
1327 for i in range(len(self)):
1328 self.__items[i]-=other[i]
1329 return self
1330
1331 def __neg__(self):
1332 """
1333 negate
1334
1335 """
1336 for i in range(len(self)):
1337 self.__items[i]=-self.__items[i]
1338 return self
1339
1340
1341 class HomogeneousSaddlePointProblem(object):
1342 """
1343 This provides a framwork for solving homogeneous saddle point problem of the form
1344
1345 Av+B^*p=f
1346 Bv =0
1347
1348 for the unknowns v and p and given operators A and B and given right hand side f.
1349 B^* is the adjoint operator of B is the given inner product.
1350
1351 """
1352 def __init__(self,**kwargs):
1353 self.setTolerance()
1354 self.setToleranceReductionFactor()
1355
1356 def initialize(self):
1357 """
1358 initialize the problem (overwrite)
1359 """
1360 pass
1361 def B(self,v):
1362 """
1363 returns Bv (overwrite)
1364 @rtype: equal to the type of p
1365
1366 @note: boundary conditions on p should be zero!
1367 """
1368 pass
1369
1370 def inner(self,p0,p1):
1371 """
1372 returns inner product of two element p0 and p1 (overwrite)
1373
1374 @type p0: equal to the type of p
1375 @type p1: equal to the type of p
1376 @rtype: C{float}
1377
1378 @rtype: equal to the type of p
1379 """
1380 pass
1381
1382 def solve_A(self,u,p):
1383 """
1384 solves Av=f-Au-B^*p with accuracy self.getReducedTolerance() (overwrite)
1385
1386 @rtype: equal to the type of v
1387 @note: boundary conditions on v should be zero!
1388 """
1389 pass
1390
1391 def solve_prec(self,p):
1392 """
1393 provides a preconditioner for BA^{-1}B^* with accuracy self.getReducedTolerance() (overwrite)
1394
1395 @rtype: equal to the type of p
1396 """
1397 pass
1398
1399 def stoppingcriterium(self,Bv,v,p):
1400 """
1401 returns a True if iteration is terminated. (overwrite)
1402
1403 @rtype: C{bool}
1404 """
1405 pass
1406
1407 def __inner(self,p,r):
1408 return self.inner(p,r[1])
1409
1410 def __inner_p(self,p1,p2):
1411 return self.inner(p1,p2)
1412
1413 def __inner_a(self,a1,a2):
1414 return self.inner_a(a1,a2)
1415
1416 def __inner_a1(self,a1,a2):
1417 return self.inner(a1[1],a2[1])
1418
1419 def __stoppingcriterium(self,norm_r,r,p):
1420 return self.stoppingcriterium(r[1],r[0],p)
1421
1422 def __stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None):
1423 return self.stoppingcriterium2(norm_r,norm_b,solver,TOL)
1424
1425 def setTolerance(self,tolerance=1.e-8):
1426 self.__tol=tolerance
1427 def getTolerance(self):
1428 return self.__tol
1429 def setToleranceReductionFactor(self,reduction=0.01):
1430 self.__reduction=reduction
1431 def getSubProblemTolerance(self):
1432 return self.__reduction*self.getTolerance()
1433
1434 def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='PCG',iter_restart=20):
1435 """
1436 solves the saddle point problem using initial guesses v and p.
1437
1438 @param max_iter: maximum number of iteration steps.
1439 """
1440 self.verbose=verbose
1441 self.show_details=show_details and self.verbose
1442
1443 # assume p is known: then v=A^-1(f-B^*p)
1444 # which leads to BA^-1B^*p = BA^-1f
1445
1446 # Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask)
1447 self.__z=v+self.solve_A(v,p*0)
1448 Bz=self.B(self.__z)
1449 #
1450 # solve BA^-1B^*p = Bz
1451 #
1452 #
1453 #
1454 self.iter=0
1455 if solver=='GMRES':
1456 if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter
1457 p=GMRES(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.,iter_restart=iter_restart)
1458 # solve Au=f-B^*p
1459 # A(u-v)=f-B^*p-Av
1460 # u=v+(u-v)
1461 u=v+self.solve_A(v,p)
1462
1463 if solver=='NewtonGMRES':
1464 if self.verbose: print "enter NewtonGMRES method (iter_max=%s)"%max_iter
1465 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())
1466 # solve Au=f-B^*p
1467 # A(u-v)=f-B^*p-Av
1468 # u=v+(u-v)
1469 u=v+self.solve_A(v,p)
1470
1471
1472 if solver=='TFQMR':
1473 if self.verbose: print "enter TFQMR method (iter_max=%s)"%max_iter
1474 p=TFQMR(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)
1475 # solve Au=f-B^*p
1476 # A(u-v)=f-B^*p-Av
1477 # u=v+(u-v)
1478 u=v+self.solve_A(v,p)
1479
1480 if solver=='MINRES':
1481 if self.verbose: print "enter MINRES method (iter_max=%s)"%max_iter
1482 p=MINRES(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)
1483 # solve Au=f-B^*p
1484 # A(u-v)=f-B^*p-Av
1485 # u=v+(u-v)
1486 u=v+self.solve_A(v,p)
1487
1488 if solver=='GMRESC':
1489 if self.verbose: print "enter GMRES coupled method (iter_max=%s)"%max_iter
1490 p0=self.solve_prec1(Bz)
1491 #alfa=(1/self.vol)*util.integrate(util.interpolate(p,escript.Function(self.domain)))
1492 #p-=alfa
1493 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)
1494 #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())
1495
1496 # solve Au=f-B^*p
1497 # A(u-v)=f-B^*p-Av
1498 # u=v+(u-v)
1499 p=x[1]
1500 u=v+self.solve_A(v,p)
1501 #p=x[1]
1502 #u=x[0]
1503
1504 if solver=='PCG':
1505 # note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv
1506 #
1507 # with Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask)
1508 # A(v-z)= f -Az - B^*p (v-z=0 on fixed_u_mask)
1509 if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter
1510 p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p)
1511 u=r[0]
1512 # print "DIFF=",util.integrate(p)
1513
1514 # print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)
1515
1516 return u,p
1517
1518 def __Msolve(self,r):
1519 return self.solve_prec(r[1])
1520
1521 def __Msolve2(self,r):
1522 return self.solve_prec(r*1)
1523
1524 def __Mempty(self,r):
1525 return r
1526
1527
1528 def __Aprod(self,p):
1529 # return BA^-1B*p
1530 #solve Av =B^*p as Av =f-Az-B^*(-p)
1531 v=self.solve_A(self.__z,-p)
1532 return ArithmeticTuple(v, self.B(v))
1533
1534 def __Anext(self,x):
1535 # return next v,p
1536 #solve Av =-B^*p as Av =f-Az-B^*p
1537
1538 pc=x[1]
1539 v=self.solve_A(self.__z,-pc)
1540 p=self.solve_prec1(self.B(v))
1541
1542 return ArithmeticTuple(v,p)
1543
1544
1545 def __Aprod2(self,p):
1546 # return BA^-1B*p
1547 #solve Av =B^*p as Av =f-Az-B^*(-p)
1548 v=self.solve_A(self.__z,-p)
1549 return self.B(v)
1550
1551 def __Aprod_Newton(self,p):
1552 # return BA^-1B*p - Bz
1553 #solve Av =-B^*p as Av =f-Az-B^*p
1554 v=self.solve_A(self.__z,-p)
1555 return self.B(v-self.__z)
1556
1557 def __Aprod_Newton2(self,x):
1558 # return BA^-1B*p - Bz
1559 #solve Av =-B^*p as Av =f-Az-B^*p
1560 pc=x[1]
1561 v=self.solve_A(self.__z,-pc)
1562 p=self.solve_prec1(self.B(v-self.__z))
1563 return ArithmeticTuple(v,p)
1564
1565 class SaddlePointProblem(object):
1566 """
1567 This implements a solver for a saddlepoint problem
1568
1569 M{f(u,p)=0}
1570 M{g(u)=0}
1571
1572 for u and p. The problem is solved with an inexact Uszawa scheme for p:
1573
1574 M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})}
1575 M{Q_g (p^{k+1}-p^{k}) = g(u^{k+1})}
1576
1577 where Q_f is an approximation of the Jacobiean A_f of f with respect to u and Q_f is an approximation of
1578 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'
1579 Q_g can be difficult, non-linear conjugate gradient method is applied to solve for p, so Q_g plays
1580 in fact the role of a preconditioner.
1581 """
1582 def __init__(self,verbose=False,*args):
1583 """
1584 initializes the problem
1585
1586 @param verbose: switches on the printing out some information
1587 @type verbose: C{bool}
1588 @note: this method may be overwritten by a particular saddle point problem
1589 """
1590 print "SaddlePointProblem should not be used anymore!"
1591 if not isinstance(verbose,bool):
1592 raise TypeError("verbose needs to be of type bool.")
1593 self.__verbose=verbose
1594 self.relaxation=1.
1595
1596 def trace(self,text):
1597 """
1598 prints text if verbose has been set
1599
1600 @param text: a text message
1601 @type text: C{str}
1602 """
1603 if self.__verbose: print "%s: %s"%(str(self),text)
1604
1605 def solve_f(self,u,p,tol=1.e-8):
1606 """
1607 solves
1608
1609 A_f du = f(u,p)
1610
1611 with tolerance C{tol} and return du. A_f is Jacobiean of f with respect to u.
1612
1613 @param u: current approximation of u
1614 @type u: L{escript.Data}
1615 @param p: current approximation of p
1616 @type p: L{escript.Data}
1617 @param tol: tolerance expected for du
1618 @type tol: C{float}
1619 @return: increment du
1620 @rtype: L{escript.Data}
1621 @note: this method has to be overwritten by a particular saddle point problem
1622 """
1623 pass
1624
1625 def solve_g(self,u,tol=1.e-8):
1626 """
1627 solves
1628
1629 Q_g dp = g(u)
1630
1631 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.
1632
1633 @param u: current approximation of u
1634 @type u: L{escript.Data}
1635 @param tol: tolerance expected for dp
1636 @type tol: C{float}
1637 @return: increment dp
1638 @rtype: L{escript.Data}
1639 @note: this method has to be overwritten by a particular saddle point problem
1640 """
1641 pass
1642
1643 def inner(self,p0,p1):
1644 """
1645 inner product of p0 and p1 approximating p. Typically this returns integrate(p0*p1)
1646 @return: inner product of p0 and p1
1647 @rtype: C{float}
1648 """
1649 pass
1650
1651 subiter_max=3
1652 def solve(self,u0,p0,tolerance=1.e-6,tolerance_u=None,iter_max=100,accepted_reduction=0.995,relaxation=None):
1653 """
1654 runs the solver
1655
1656 @param u0: initial guess for C{u}
1657 @type u0: L{esys.escript.Data}
1658 @param p0: initial guess for C{p}
1659 @type p0: L{esys.escript.Data}
1660 @param tolerance: tolerance for relative error in C{u} and C{p}
1661 @type tolerance: positive C{float}
1662 @param tolerance_u: tolerance for relative error in C{u} if different from C{tolerance}
1663 @type tolerance_u: positive C{float}
1664 @param iter_max: maximum number of iteration steps.
1665 @type iter_max: C{int}
1666 @param accepted_reduction: if the norm g cannot be reduced by C{accepted_reduction} backtracking to adapt the
1667 relaxation factor. If C{accepted_reduction=None} no backtracking is used.
1668 @type accepted_reduction: positive C{float} or C{None}
1669 @param relaxation: initial relaxation factor. If C{relaxation==None}, the last relaxation factor is used.
1670 @type relaxation: C{float} or C{None}
1671 """
1672 tol=1.e-2
1673 if tolerance_u==None: tolerance_u=tolerance
1674 if not relaxation==None: self.relaxation=relaxation
1675 if accepted_reduction ==None:
1676 angle_limit=0.
1677 elif accepted_reduction>=1.:
1678 angle_limit=0.
1679 else:
1680 angle_limit=util.sqrt(1-accepted_reduction**2)
1681 self.iter=0
1682 u=u0
1683 p=p0
1684 #
1685 # initialize things:
1686 #
1687 converged=False
1688 #
1689 # start loop:
1690 #
1691 # initial search direction is g
1692 #
1693 while not converged :
1694 if self.iter>iter_max:
1695 raise ArithmeticError("no convergence after %s steps."%self.iter)
1696 f_new=self.solve_f(u,p,tol)
1697 norm_f_new = util.Lsup(f_new)
1698 u_new=u-f_new
1699 g_new=self.solve_g(u_new,tol)
1700 self.iter+=1
1701 norm_g_new = util.sqrt(self.inner(g_new,g_new))
1702 if norm_f_new==0. and norm_g_new==0.: return u, p
1703 if self.iter>1 and not accepted_reduction==None:
1704 #
1705 # did we manage to reduce the norm of G? I
1706 # if not we start a backtracking procedure
1707 #
1708 # print "new/old norm = ",norm_g_new, norm_g, norm_g_new/norm_g
1709 if norm_g_new > accepted_reduction * norm_g:
1710 sub_iter=0
1711 s=self.relaxation
1712 d=g
1713 g_last=g
1714 self.trace(" start substepping: f = %s, g = %s, relaxation = %s."%(norm_f_new, norm_g_new, s))
1715 while sub_iter < self.subiter_max and norm_g_new > accepted_reduction * norm_g:
1716 dg= g_new-g_last
1717 norm_dg=abs(util.sqrt(self.inner(dg,dg))/self.relaxation)
1718 rad=self.inner(g_new,dg)/self.relaxation
1719 # print " ",sub_iter,": rad, norm_dg:",abs(rad), norm_dg*norm_g_new * angle_limit
1720 # print " ",sub_iter,": rad, norm_dg:",rad, norm_dg, norm_g_new, norm_g
1721 if abs(rad) < norm_dg*norm_g_new * angle_limit:
1722 if sub_iter>0: self.trace(" no further improvements expected from backtracking.")
1723 break
1724 r=self.relaxation
1725 self.relaxation= - rad/norm_dg**2
1726 s+=self.relaxation
1727 #####
1728 # a=g_new+self.relaxation*dg/r
1729 # print "predicted new norm = ",util.sqrt(self.inner(a,a)),util.sqrt(self.inner(g_new,g_new)), self.relaxation
1730 #####
1731 g_last=g_new
1732 p+=self.relaxation*d
1733 f_new=self.solve_f(u,p,tol)
1734 u_new=u-f_new
1735 g_new=self.solve_g(u_new,tol)
1736 self.iter+=1
1737 norm_f_new = util.Lsup(f_new)
1738 norm_g_new = util.sqrt(self.inner(g_new,g_new))
1739 # print " ",sub_iter," new g norm",norm_g_new
1740 self.trace(" %s th sub-step: f = %s, g = %s, relaxation = %s."%(sub_iter, norm_f_new, norm_g_new, s))
1741 #
1742 # can we expect reduction of g?
1743 #
1744 # u_last=u_new
1745 sub_iter+=1
1746 self.relaxation=s
1747 #
1748 # check for convergence:
1749 #
1750 norm_u_new = util.Lsup(u_new)
1751 p_new=p+self.relaxation*g_new
1752 norm_p_new = util.sqrt(self.inner(p_new,p_new))
1753 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))
1754
1755 if self.iter>1:
1756 dg2=g_new-g
1757 df2=f_new-f
1758 norm_dg2=util.sqrt(self.inner(dg2,dg2))
1759 norm_df2=util.Lsup(df2)
1760 # print norm_g_new, norm_g, norm_dg, norm_p, tolerance
1761 tol_eq_g=tolerance*norm_dg2/(norm_g*abs(self.relaxation))*norm_p_new
1762 tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new
1763 if norm_g_new <= tol_eq_g and norm_f_new <= tol_eq_f:
1764 converged=True
1765 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
1766 self.trace("convergence after %s steps."%self.iter)
1767 return u,p
1768 # def solve(self,u0,p0,tolerance=1.e-6,iter_max=10,self.relaxation=1.):
1769 # tol=1.e-2
1770 # iter=0
1771 # converged=False
1772 # u=u0*1.
1773 # p=p0*1.
1774 # while not converged and iter<iter_max:
1775 # du=self.solve_f(u,p,tol)
1776 # u-=du
1777 # norm_du=util.Lsup(du)
1778 # norm_u=util.Lsup(u)
1779 #
1780 # dp=self.relaxation*self.solve_g(u,tol)
1781 # p+=dp
1782 # norm_dp=util.sqrt(self.inner(dp,dp))
1783 # norm_p=util.sqrt(self.inner(p,p))
1784 # 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)
1785 # iter+=1
1786 #
1787 # converged = (norm_du <= tolerance*norm_u) and (norm_dp <= tolerance*norm_p)
1788 # if converged:
1789 # print "convergence after %s steps."%iter
1790 # else:
1791 # raise ArithmeticError("no convergence after %s steps."%iter)
1792 #
1793 # return u,p
1794
1795 def MaskFromBoundaryTag(domain,*tags):
1796 """
1797 create a mask on the Solution(domain) function space which one for samples
1798 that touch the boundary tagged by tags.
1799
1800 usage: m=MaskFromBoundaryTag(domain,"left", "right")
1801
1802 @param domain: a given domain
1803 @type domain: L{escript.Domain}
1804 @param tags: boundray tags
1805 @type tags: C{str}
1806 @return: a mask which marks samples that are touching the boundary tagged by any of the given tags.
1807 @rtype: L{escript.Data} of rank 0
1808 """
1809 pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1810 d=escript.Scalar(0.,escript.FunctionOnBoundary(function_space.getDomain()))
1811 for t in tags: d.setTaggedValue(t,1.)
1812 pde.setValue(y=d)
1813 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