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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26