/[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 1557 - (show annotations)
Mon May 19 04:44:27 2008 UTC (11 years, 4 months ago) by artak
File MIME type: text/x-python
File size: 57307 byte(s)
some changes to ArithmenticTuple class
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
1048 etamax=min(etanew,etamax)
1049 etamax=max(etamax,.5*stop_tol/fnrm)
1050
1051 return x
1052
1053 def TFQMR(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
1054
1055 # TFQMR solver for linear systems
1056 #
1057 #
1058 # initialization
1059 #
1060 errtol = math.sqrt(bilinearform(b,b))
1061 norm_b=errtol
1062 kmax = iter_max
1063 error = []
1064
1065 if math.sqrt(bilinearform(x,x)) != 0.0:
1066 r = b - Aprod(x)
1067 else:
1068 r = b
1069
1070 r=Msolve(r)
1071
1072 u1=0
1073 u2=0
1074 y1=0
1075 y2=0
1076
1077 w = r
1078 y1 = r
1079 iter = 0
1080 d = 0
1081
1082 v = Msolve(Aprod(y1))
1083 u1 = v
1084
1085 theta = 0.0;
1086 eta = 0.0;
1087 tau = math.sqrt(bilinearform(r,r))
1088 error = [ error, tau ]
1089 rho = tau * tau
1090 m=1
1091 #
1092 # TFQMR iteration
1093 #
1094 # while ( iter < kmax-1 ):
1095
1096 while not stoppingcriterium(tau*math.sqrt ( m + 1 ),norm_b,'TFQMR'):
1097 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
1098
1099 sigma = bilinearform(r,v)
1100
1101 if ( sigma == 0.0 ):
1102 raise 'TFQMR breakdown, sigma=0'
1103
1104
1105 alpha = rho / sigma
1106
1107 for j in range(2):
1108 #
1109 # Compute y2 and u2 only if you have to
1110 #
1111 if ( j == 1 ):
1112 y2 = y1 - alpha * v
1113 u2 = Msolve(Aprod(y2))
1114
1115 m = 2 * (iter+1) - 2 + (j+1)
1116 if j==0:
1117 w = w - alpha * u1
1118 d = y1 + ( theta * theta * eta / alpha ) * d
1119 if j==1:
1120 w = w - alpha * u2
1121 d = y2 + ( theta * theta * eta / alpha ) * d
1122
1123 theta = math.sqrt(bilinearform(w,w))/ tau
1124 c = 1.0 / math.sqrt ( 1.0 + theta * theta )
1125 tau = tau * theta * c
1126 eta = c * c * alpha
1127 x = x + eta * d
1128 #
1129 # Try to terminate the iteration at each pass through the loop
1130 #
1131 # if ( tau * math.sqrt ( m + 1 ) <= errtol ):
1132 # error = [ error, tau ]
1133 # total_iters = iter
1134 # break
1135
1136
1137 if ( rho == 0.0 ):
1138 raise 'TFQMR breakdown, rho=0'
1139
1140
1141 rhon = bilinearform(r,w)
1142 beta = rhon / rho;
1143 rho = rhon;
1144 y1 = w + beta * y2;
1145 u1 = Msolve(Aprod(y1))
1146 v = u1 + beta * ( u2 + beta * v )
1147 error = [ error, tau ]
1148 total_iters = iter
1149
1150 iter = iter + 1
1151
1152 return x
1153
1154
1155 #############################################
1156
1157 class ArithmeticTuple(object):
1158 """
1159 tuple supporting inplace update x+=y and scaling x=a*y where x,y is an ArithmeticTuple and a is a float.
1160
1161 example of usage:
1162
1163 from esys.escript import Data
1164 from numarray import array
1165 a=Data(...)
1166 b=array([1.,4.])
1167 x=ArithmeticTuple(a,b)
1168 y=5.*x
1169
1170 """
1171 def __init__(self,*args):
1172 """
1173 initialize object with elements args.
1174
1175 @param args: tuple of object that support implace add (x+=y) and scaling (x=a*y)
1176 """
1177 self.__items=list(args)
1178
1179 def __len__(self):
1180 """
1181 number of items
1182
1183 @return: number of items
1184 @rtype: C{int}
1185 """
1186 return len(self.__items)
1187
1188 def __getitem__(self,index):
1189 """
1190 get an item
1191
1192 @param index: item to be returned
1193 @type index: C{int}
1194 @return: item with index C{index}
1195 """
1196 return self.__items.__getitem__(index)
1197
1198 def __mul__(self,other):
1199 """
1200 scaling from the right
1201
1202 @param other: scaling factor
1203 @type other: C{float}
1204 @return: itemwise self*other
1205 @rtype: L{ArithmeticTuple}
1206 """
1207 out=[]
1208 other=1.*other
1209 if isinstance(other,float):
1210 for i in range(len(self)):
1211 out.append(self[i]*other)
1212 else:
1213 for i in range(len(self)):
1214 out.append(self[i]*other[i])
1215 return ArithmeticTuple(*tuple(out))
1216
1217 def __rmul__(self,other):
1218 """
1219 scaling from the left
1220
1221 @param other: scaling factor
1222 @type other: C{float}
1223 @return: itemwise other*self
1224 @rtype: L{ArithmeticTuple}
1225 """
1226 out=[]
1227 other=1.*other
1228 if isinstance(other,float):
1229 for i in range(len(self)):
1230 out.append(other*self[i])
1231 else:
1232 for i in range(len(self)):
1233 out.append(other[i]*self[i])
1234 return ArithmeticTuple(*tuple(out))
1235
1236 #########################
1237 # Added by Artak
1238 #########################
1239 def __div__(self,other):
1240 """
1241 dividing from the right
1242
1243 @param other: scaling factor
1244 @type other: C{float}
1245 @return: itemwise self/other
1246 @rtype: L{ArithmeticTuple}
1247 """
1248 out=[]
1249 other=1.*other
1250 if isinstance(other,float):
1251 for i in range(len(self)):
1252 out.append(self[i]*(1./other))
1253 else:
1254 for i in range(len(self)):
1255 out.append(self[i]*(1./other[i]))
1256 return ArithmeticTuple(*tuple(out))
1257
1258 def __rdiv__(self,other):
1259 """
1260 dividing from the left
1261
1262 @param other: scaling factor
1263 @type other: C{float}
1264 @return: itemwise other/self
1265 @rtype: L{ArithmeticTuple}
1266 """
1267 out=[]
1268 other=1.*other
1269 if isinstance(other,float):
1270 for i in range(len(self)):
1271 out.append(other*(1./self[i]))
1272 else:
1273 for i in range(len(self)):
1274 out.append(other[i]*(1./self[i]))
1275 return ArithmeticTuple(*tuple(out))
1276
1277 ##########################################33
1278
1279 def __iadd__(self,other):
1280 """
1281 in-place add of other to self
1282
1283 @param other: increment
1284 @type other: C{ArithmeticTuple}
1285 """
1286 if len(self) != len(other):
1287 raise ValueError,"tuple length must match."
1288 for i in range(len(self)):
1289 self.__items[i]+=other[i]
1290 return self
1291
1292 def __add__(self,other):
1293 """
1294 add other to self
1295
1296 @param other: increment
1297 @type other: C{ArithmeticTuple}
1298 """
1299 # if not isinstance(other,float):
1300 if len(self) != len(other):
1301 raise ValueError,"tuple length must match."
1302 for i in range(len(self)):
1303 self.__items[i]+=other[i]
1304 # else:
1305 # for i in range(len(self)):
1306 # self.__items[i]+=other
1307
1308 return self
1309
1310 def __sub__(self,other):
1311 """
1312 subtract other from self
1313
1314 @param other: increment
1315 @type other: C{ArithmeticTuple}
1316 """
1317 if len(self) != len(other):
1318 raise ValueError,"tuple length must match."
1319 for i in range(len(self)):
1320 self.__items[i]-=other[i]
1321 return self
1322
1323 def __isub__(self,other):
1324 """
1325 subtract other from self
1326
1327 @param other: increment
1328 @type other: C{ArithmeticTuple}
1329 """
1330 if len(self) != len(other):
1331 raise ValueError,"tuple length must match."
1332 for i in range(len(self)):
1333 self.__items[i]-=other[i]
1334 return self
1335
1336 def __neg__(self):
1337 """
1338 negate
1339
1340 """
1341 for i in range(len(self)):
1342 self.__items[i]=-self.__items[i]
1343 return self
1344
1345
1346 class HomogeneousSaddlePointProblem(object):
1347 """
1348 This provides a framwork for solving homogeneous saddle point problem of the form
1349
1350 Av+B^*p=f
1351 Bv =0
1352
1353 for the unknowns v and p and given operators A and B and given right hand side f.
1354 B^* is the adjoint operator of B is the given inner product.
1355
1356 """
1357 def __init__(self,**kwargs):
1358 self.setTolerance()
1359 self.setToleranceReductionFactor()
1360
1361 def initialize(self):
1362 """
1363 initialize the problem (overwrite)
1364 """
1365 pass
1366 def B(self,v):
1367 """
1368 returns Bv (overwrite)
1369 @rtype: equal to the type of p
1370
1371 @note: boundary conditions on p should be zero!
1372 """
1373 pass
1374
1375 def inner(self,p0,p1):
1376 """
1377 returns inner product of two element p0 and p1 (overwrite)
1378
1379 @type p0: equal to the type of p
1380 @type p1: equal to the type of p
1381 @rtype: C{float}
1382
1383 @rtype: equal to the type of p
1384 """
1385 pass
1386
1387 def solve_A(self,u,p):
1388 """
1389 solves Av=f-Au-B^*p with accuracy self.getReducedTolerance() (overwrite)
1390
1391 @rtype: equal to the type of v
1392 @note: boundary conditions on v should be zero!
1393 """
1394 pass
1395
1396 def solve_prec(self,p):
1397 """
1398 provides a preconditioner for BA^{-1}B^* with accuracy self.getReducedTolerance() (overwrite)
1399
1400 @rtype: equal to the type of p
1401 """
1402 pass
1403
1404 def stoppingcriterium(self,Bv,v,p):
1405 """
1406 returns a True if iteration is terminated. (overwrite)
1407
1408 @rtype: C{bool}
1409 """
1410 pass
1411
1412 def __inner(self,p,r):
1413 return self.inner(p,r[1])
1414
1415 def __inner_p(self,p1,p2):
1416 return self.inner(p1,p2)
1417
1418 def __inner_a(self,a1,a2):
1419 return self.inner_a(a1,a2)
1420
1421 def __inner_a1(self,a1,a2):
1422 return self.inner(a1[1],a2[1])
1423
1424 def __stoppingcriterium(self,norm_r,r,p):
1425 return self.stoppingcriterium(r[1],r[0],p)
1426
1427 def __stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None):
1428 return self.stoppingcriterium2(norm_r,norm_b,solver,TOL)
1429
1430 def setTolerance(self,tolerance=1.e-8):
1431 self.__tol=tolerance
1432 def getTolerance(self):
1433 return self.__tol
1434 def setToleranceReductionFactor(self,reduction=0.01):
1435 self.__reduction=reduction
1436 def getSubProblemTolerance(self):
1437 return self.__reduction*self.getTolerance()
1438
1439 def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='PCG',iter_restart=20):
1440 """
1441 solves the saddle point problem using initial guesses v and p.
1442
1443 @param max_iter: maximum number of iteration steps.
1444 """
1445 self.verbose=verbose
1446 self.show_details=show_details and self.verbose
1447
1448 # assume p is known: then v=A^-1(f-B^*p)
1449 # which leads to BA^-1B^*p = BA^-1f
1450
1451 # Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask)
1452 self.__z=v+self.solve_A(v,p*0)
1453 Bz=self.B(self.__z)
1454 #
1455 # solve BA^-1B^*p = Bz
1456 #
1457 #
1458 #
1459 self.iter=0
1460 if solver=='GMRES':
1461 if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter
1462 p=GMRES(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.,iter_restart=iter_restart)
1463 # solve Au=f-B^*p
1464 # A(u-v)=f-B^*p-Av
1465 # u=v+(u-v)
1466 u=v+self.solve_A(v,p)
1467
1468 if solver=='NewtonGMRES':
1469 if self.verbose: print "enter NewtonGMRES method (iter_max=%s)"%max_iter
1470 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())
1471 # solve Au=f-B^*p
1472 # A(u-v)=f-B^*p-Av
1473 # u=v+(u-v)
1474 u=v+self.solve_A(v,p)
1475
1476
1477 if solver=='TFQMR':
1478 if self.verbose: print "enter TFQMR method (iter_max=%s)"%max_iter
1479 p=TFQMR(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)
1480 # solve Au=f-B^*p
1481 # A(u-v)=f-B^*p-Av
1482 # u=v+(u-v)
1483 u=v+self.solve_A(v,p)
1484
1485 if solver=='MINRES':
1486 if self.verbose: print "enter MINRES method (iter_max=%s)"%max_iter
1487 p=MINRES(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)
1488 # solve Au=f-B^*p
1489 # A(u-v)=f-B^*p-Av
1490 # u=v+(u-v)
1491 u=v+self.solve_A(v,p)
1492
1493 if solver=='GMRESC':
1494 if self.verbose: print "enter GMRES coupled method (iter_max=%s)"%max_iter
1495 p0=self.solve_prec1(Bz)
1496 #alfa=(1/self.vol)*util.integrate(util.interpolate(p,escript.Function(self.domain)))
1497 #p-=alfa
1498 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)
1499 #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())
1500
1501 # solve Au=f-B^*p
1502 # A(u-v)=f-B^*p-Av
1503 # u=v+(u-v)
1504 p=x[1]
1505 u=v+self.solve_A(v,p)
1506 #p=x[1]
1507 #u=x[0]
1508
1509 if solver=='PCG':
1510 # note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv
1511 #
1512 # with Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask)
1513 # A(v-z)= f -Az - B^*p (v-z=0 on fixed_u_mask)
1514 if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter
1515 p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p)
1516 u=r[0]
1517 print "DIFF=",util.integrate(p)
1518
1519 print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)
1520
1521 return u,p
1522
1523 def __Msolve(self,r):
1524 return self.solve_prec(r[1])
1525
1526 def __Msolve2(self,r):
1527 return self.solve_prec(r*1)
1528
1529 def __Mempty(self,r):
1530 return r
1531
1532
1533 def __Aprod(self,p):
1534 # return BA^-1B*p
1535 #solve Av =B^*p as Av =f-Az-B^*(-p)
1536 v=self.solve_A(self.__z,-p)
1537 return ArithmeticTuple(v, self.B(v))
1538
1539 def __Anext(self,x):
1540 # return next v,p
1541 #solve Av =-B^*p as Av =f-Az-B^*p
1542
1543 pc=x[1]
1544 v=self.solve_A(self.__z,-pc)
1545 p=self.solve_prec1(self.B(v))
1546
1547 return ArithmeticTuple(v,p)
1548
1549
1550 def __Aprod2(self,p):
1551 # return BA^-1B*p
1552 #solve Av =B^*p as Av =f-Az-B^*(-p)
1553 v=self.solve_A(self.__z,-p)
1554 return self.B(v)
1555
1556 def __Aprod_Newton(self,p):
1557 # return BA^-1B*p - Bz
1558 #solve Av =-B^*p as Av =f-Az-B^*p
1559 v=self.solve_A(self.__z,-p)
1560 return self.B(v-self.__z)
1561
1562 def __Aprod_Newton2(self,x):
1563 # return BA^-1B*p - Bz
1564 #solve Av =-B^*p as Av =f-Az-B^*p
1565 pc=x[1]
1566 v=self.solve_A(self.__z,-pc)
1567 p=self.solve_prec1(self.B(v-self.__z))
1568 return ArithmeticTuple(v,p)
1569
1570 class SaddlePointProblem(object):
1571 """
1572 This implements a solver for a saddlepoint problem
1573
1574 M{f(u,p)=0}
1575 M{g(u)=0}
1576
1577 for u and p. The problem is solved with an inexact Uszawa scheme for p:
1578
1579 M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})}
1580 M{Q_g (p^{k+1}-p^{k}) = g(u^{k+1})}
1581
1582 where Q_f is an approximation of the Jacobiean A_f of f with respect to u and Q_f is an approximation of
1583 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'
1584 Q_g can be difficult, non-linear conjugate gradient method is applied to solve for p, so Q_g plays
1585 in fact the role of a preconditioner.
1586 """
1587 def __init__(self,verbose=False,*args):
1588 """
1589 initializes the problem
1590
1591 @param verbose: switches on the printing out some information
1592 @type verbose: C{bool}
1593 @note: this method may be overwritten by a particular saddle point problem
1594 """
1595 if not isinstance(verbose,bool):
1596 raise TypeError("verbose needs to be of type bool.")
1597 self.__verbose=verbose
1598 self.relaxation=1.
1599
1600 def trace(self,text):
1601 """
1602 prints text if verbose has been set
1603
1604 @param text: a text message
1605 @type text: C{str}
1606 """
1607 if self.__verbose: print "#s: #s"%(str(self),text)
1608
1609 def solve_f(self,u,p,tol=1.e-8):
1610 """
1611 solves
1612
1613 A_f du = f(u,p)
1614
1615 with tolerance C{tol} and return du. A_f is Jacobiean of f with respect to u.
1616
1617 @param u: current approximation of u
1618 @type u: L{escript.Data}
1619 @param p: current approximation of p
1620 @type p: L{escript.Data}
1621 @param tol: tolerance expected for du
1622 @type tol: C{float}
1623 @return: increment du
1624 @rtype: L{escript.Data}
1625 @note: this method has to be overwritten by a particular saddle point problem
1626 """
1627 pass
1628
1629 def solve_g(self,u,tol=1.e-8):
1630 """
1631 solves
1632
1633 Q_g dp = g(u)
1634
1635 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.
1636
1637 @param u: current approximation of u
1638 @type u: L{escript.Data}
1639 @param tol: tolerance expected for dp
1640 @type tol: C{float}
1641 @return: increment dp
1642 @rtype: L{escript.Data}
1643 @note: this method has to be overwritten by a particular saddle point problem
1644 """
1645 pass
1646
1647 def inner(self,p0,p1):
1648 """
1649 inner product of p0 and p1 approximating p. Typically this returns integrate(p0*p1)
1650 @return: inner product of p0 and p1
1651 @rtype: C{float}
1652 """
1653 pass
1654
1655 subiter_max=3
1656 def solve(self,u0,p0,tolerance=1.e-6,tolerance_u=None,iter_max=100,accepted_reduction=0.995,relaxation=None):
1657 """
1658 runs the solver
1659
1660 @param u0: initial guess for C{u}
1661 @type u0: L{esys.escript.Data}
1662 @param p0: initial guess for C{p}
1663 @type p0: L{esys.escript.Data}
1664 @param tolerance: tolerance for relative error in C{u} and C{p}
1665 @type tolerance: positive C{float}
1666 @param tolerance_u: tolerance for relative error in C{u} if different from C{tolerance}
1667 @type tolerance_u: positive C{float}
1668 @param iter_max: maximum number of iteration steps.
1669 @type iter_max: C{int}
1670 @param accepted_reduction: if the norm g cannot be reduced by C{accepted_reduction} backtracking to adapt the
1671 relaxation factor. If C{accepted_reduction=None} no backtracking is used.
1672 @type accepted_reduction: positive C{float} or C{None}
1673 @param relaxation: initial relaxation factor. If C{relaxation==None}, the last relaxation factor is used.
1674 @type relaxation: C{float} or C{None}
1675 """
1676 tol=1.e-2
1677 if tolerance_u==None: tolerance_u=tolerance
1678 if not relaxation==None: self.relaxation=relaxation
1679 if accepted_reduction ==None:
1680 angle_limit=0.
1681 elif accepted_reduction>=1.:
1682 angle_limit=0.
1683 else:
1684 angle_limit=util.sqrt(1-accepted_reduction**2)
1685 self.iter=0
1686 u=u0
1687 p=p0
1688 #
1689 # initialize things:
1690 #
1691 converged=False
1692 #
1693 # start loop:
1694 #
1695 # initial search direction is g
1696 #
1697 while not converged :
1698 if self.iter>iter_max:
1699 raise ArithmeticError("no convergence after %s steps."%self.iter)
1700 f_new=self.solve_f(u,p,tol)
1701 norm_f_new = util.Lsup(f_new)
1702 u_new=u-f_new
1703 g_new=self.solve_g(u_new,tol)
1704 self.iter+=1
1705 norm_g_new = util.sqrt(self.inner(g_new,g_new))
1706 if norm_f_new==0. and norm_g_new==0.: return u, p
1707 if self.iter>1 and not accepted_reduction==None:
1708 #
1709 # did we manage to reduce the norm of G? I
1710 # if not we start a backtracking procedure
1711 #
1712 # print "new/old norm = ",norm_g_new, norm_g, norm_g_new/norm_g
1713 if norm_g_new > accepted_reduction * norm_g:
1714 sub_iter=0
1715 s=self.relaxation
1716 d=g
1717 g_last=g
1718 self.trace(" start substepping: f = %s, g = %s, relaxation = %s."%(norm_f_new, norm_g_new, s))
1719 while sub_iter < self.subiter_max and norm_g_new > accepted_reduction * norm_g:
1720 dg= g_new-g_last
1721 norm_dg=abs(util.sqrt(self.inner(dg,dg))/self.relaxation)
1722 rad=self.inner(g_new,dg)/self.relaxation
1723 # print " ",sub_iter,": rad, norm_dg:",abs(rad), norm_dg*norm_g_new * angle_limit
1724 # print " ",sub_iter,": rad, norm_dg:",rad, norm_dg, norm_g_new, norm_g
1725 if abs(rad) < norm_dg*norm_g_new * angle_limit:
1726 if sub_iter>0: self.trace(" no further improvements expected from backtracking.")
1727 break
1728 r=self.relaxation
1729 self.relaxation= - rad/norm_dg**2
1730 s+=self.relaxation
1731 #####
1732 # a=g_new+self.relaxation*dg/r
1733 # print "predicted new norm = ",util.sqrt(self.inner(a,a)),util.sqrt(self.inner(g_new,g_new)), self.relaxation
1734 #####
1735 g_last=g_new
1736 p+=self.relaxation*d
1737 f_new=self.solve_f(u,p,tol)
1738 u_new=u-f_new
1739 g_new=self.solve_g(u_new,tol)
1740 self.iter+=1
1741 norm_f_new = util.Lsup(f_new)
1742 norm_g_new = util.sqrt(self.inner(g_new,g_new))
1743 # print " ",sub_iter," new g norm",norm_g_new
1744 self.trace(" %s th sub-step: f = %s, g = %s, relaxation = %s."%(sub_iter, norm_f_new, norm_g_new, s))
1745 #
1746 # can we expect reduction of g?
1747 #
1748 # u_last=u_new
1749 sub_iter+=1
1750 self.relaxation=s
1751 #
1752 # check for convergence:
1753 #
1754 norm_u_new = util.Lsup(u_new)
1755 p_new=p+self.relaxation*g_new
1756 norm_p_new = util.sqrt(self.inner(p_new,p_new))
1757 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))
1758
1759 if self.iter>1:
1760 dg2=g_new-g
1761 df2=f_new-f
1762 norm_dg2=util.sqrt(self.inner(dg2,dg2))
1763 norm_df2=util.Lsup(df2)
1764 # print norm_g_new, norm_g, norm_dg, norm_p, tolerance
1765 tol_eq_g=tolerance*norm_dg2/(norm_g*abs(self.relaxation))*norm_p_new
1766 tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new
1767 if norm_g_new <= tol_eq_g and norm_f_new <= tol_eq_f:
1768 converged=True
1769 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
1770 self.trace("convergence after %s steps."%self.iter)
1771 return u,p
1772 # def solve(self,u0,p0,tolerance=1.e-6,iter_max=10,self.relaxation=1.):
1773 # tol=1.e-2
1774 # iter=0
1775 # converged=False
1776 # u=u0*1.
1777 # p=p0*1.
1778 # while not converged and iter<iter_max:
1779 # du=self.solve_f(u,p,tol)
1780 # u-=du
1781 # norm_du=util.Lsup(du)
1782 # norm_u=util.Lsup(u)
1783 #
1784 # dp=self.relaxation*self.solve_g(u,tol)
1785 # p+=dp
1786 # norm_dp=util.sqrt(self.inner(dp,dp))
1787 # norm_p=util.sqrt(self.inner(p,p))
1788 # 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)
1789 # iter+=1
1790 #
1791 # converged = (norm_du <= tolerance*norm_u) and (norm_dp <= tolerance*norm_p)
1792 # if converged:
1793 # print "convergence after %s steps."%iter
1794 # else:
1795 # raise ArithmeticError("no convergence after %s steps."%iter)
1796 #
1797 # return u,p
1798
1799 def MaskFromBoundaryTag(function_space,*tags):
1800 """
1801 create a mask on the given function space which one for samples
1802 that touch the boundary tagged by tags.
1803
1804 usage: m=MaskFromBoundaryTag(Solution(domain),"left", "right")
1805
1806 @param function_space: a given function space
1807 @type function_space: L{escript.FunctionSpace}
1808 @param tags: boundray tags
1809 @type tags: C{str}
1810 @return: a mask which marks samples used by C{function_space} that are touching the
1811 boundary tagged by any of the given tags.
1812 @rtype: L{escript.Data} of rank 0
1813 """
1814 pde=linearPDEs.LinearPDE(function_space.getDomain(),numEquations=1, numSolutions=1)
1815 d=escript.Scalar(0.,escript.FunctionOnBoundary(function_space.getDomain()))
1816 for t in tags: d.setTaggedValue(t,1.)
1817 pde.setValue(y=d)
1818 out=util.whereNonZero(pde.getRightHandSide())
1819 if out.getFunctionSpace() == function_space:
1820 return out
1821 else:
1822 return util.whereNonZero(util.interpolate(out,function_space))
1823
1824
1825

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26