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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26