/[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 2303 - (show annotations)
Fri Mar 13 01:34:36 2009 UTC (10 years, 8 months ago) by artak
File MIME type: text/x-python
File size: 64344 byte(s)
This was the troublemaker for test failures on both MAC and Debian. I think numarry has a bug in memory managment. So we removed some numarray.array related lines in _givapp function and tests pass again.
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,math.sqrt(rhat_dot_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]=w1
732 vrot[i+1]=w2
733 return vrot
734
735 def __FDGMRES(F0, defect, x0, atol, iter_max=100, iter_restart=20):
736 h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)
737 c=numarray.zeros(iter_restart,numarray.Float64)
738 s=numarray.zeros(iter_restart,numarray.Float64)
739 g=numarray.zeros(iter_restart,numarray.Float64)
740 v=[]
741
742 rho=defect.norm(F0)
743 if rho<=0.: return x0*0
744
745 v.append(-F0/rho)
746 g[0]=rho
747 iter=0
748 while rho > atol and iter<iter_restart-1:
749 if iter >= iter_max:
750 raise MaxIterReached,"maximum number of %s steps reached."%iter_max
751
752 p=defect.derivative(F0,x0,v[iter], v_is_normalised=True)
753 v.append(p)
754
755 v_norm1=defect.norm(v[iter+1])
756
757 # Modified Gram-Schmidt
758 for j in range(iter+1):
759 h[j,iter]=defect.bilinearform(v[j],v[iter+1])
760 v[iter+1]-=h[j,iter]*v[j]
761
762 h[iter+1,iter]=defect.norm(v[iter+1])
763 v_norm2=h[iter+1,iter]
764
765 # Reorthogonalize if needed
766 if v_norm1 + 0.001*v_norm2 == v_norm1: #Brown/Hindmarsh condition (default)
767 for j in range(iter+1):
768 hr=defect.bilinearform(v[j],v[iter+1])
769 h[j,iter]=h[j,iter]+hr
770 v[iter+1] -= hr*v[j]
771
772 v_norm2=defect.norm(v[iter+1])
773 h[iter+1,iter]=v_norm2
774 # watch out for happy breakdown
775 if not v_norm2 == 0:
776 v[iter+1]=v[iter+1]/h[iter+1,iter]
777
778 # Form and store the information for the new Givens rotation
779 if iter > 0 :
780 hhat=numarray.zeros(iter+1,numarray.Float64)
781 for i in range(iter+1) : hhat[i]=h[i,iter]
782 hhat=__givapp(c[0:iter],s[0:iter],hhat);
783 for i in range(iter+1) : h[i,iter]=hhat[i]
784
785 mu=math.sqrt(h[iter,iter]*h[iter,iter]+h[iter+1,iter]*h[iter+1,iter])
786
787 if mu!=0 :
788 c[iter]=h[iter,iter]/mu
789 s[iter]=-h[iter+1,iter]/mu
790 h[iter,iter]=c[iter]*h[iter,iter]-s[iter]*h[iter+1,iter]
791 h[iter+1,iter]=0.0
792 gg=__givapp(c[iter],s[iter],[g[iter],g[iter+1]])
793 g[iter]=gg[0]
794 g[iter+1]=gg[1]
795
796 # Update the residual norm
797 rho=abs(g[iter+1])
798 iter+=1
799
800 # At this point either iter > iter_max or rho < tol.
801 # It's time to compute x and leave.
802 if iter > 0 :
803 y=numarray.zeros(iter,numarray.Float64)
804 y[iter-1] = g[iter-1] / h[iter-1,iter-1]
805 if iter > 1 :
806 i=iter-2
807 while i>=0 :
808 y[i] = ( g[i] - numarray.dot(h[i,i+1:iter], y[i+1:iter])) / h[i,i]
809 i=i-1
810 xhat=v[iter-1]*y[iter-1]
811 for i in range(iter-1):
812 xhat += v[i]*y[i]
813 else :
814 xhat=v[0] * 0
815
816 if iter<iter_restart-1:
817 stopped=iter
818 else:
819 stopped=-1
820
821 return xhat,stopped
822
823 def GMRES(r, Aprod, x, bilinearform, atol=0, rtol=1.e-8, iter_max=100, iter_restart=20, verbose=False,P_R=None):
824 """
825 Solver for
826
827 M{Ax=b}
828
829 with a general operator A (more details required!).
830 It uses the generalized minimum residual method (GMRES).
831
832 The iteration is terminated if
833
834 M{|r| <= atol+rtol*|r0|}
835
836 where M{r0} is the initial residual and M{|.|} is the energy norm. In fact
837
838 M{|r| = sqrt( bilinearform(r,r))}
839
840 @param r: initial residual M{r=b-Ax}. C{r} is altered.
841 @type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
842 @param x: an initial guess for the solution
843 @type x: same like C{r}
844 @param Aprod: returns the value Ax
845 @type Aprod: function C{Aprod(x)} where C{x} is of the same object like
846 argument C{x}. The returned object needs to be of the same
847 type like argument C{r}.
848 @param bilinearform: inner product C{<x,r>}
849 @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same
850 type like argument C{x} and C{r}. The returned value is
851 a C{float}.
852 @param atol: absolute tolerance
853 @type atol: non-negative C{float}
854 @param rtol: relative tolerance
855 @type rtol: non-negative C{float}
856 @param iter_max: maximum number of iteration steps
857 @type iter_max: C{int}
858 @param iter_restart: in order to save memory the orthogonalization process
859 is terminated after C{iter_restart} steps and the
860 iteration is restarted.
861 @type iter_restart: C{int}
862 @return: the solution approximation and the corresponding residual
863 @rtype: C{tuple}
864 @warning: C{r} and C{x} are altered.
865 """
866 m=iter_restart
867 restarted=False
868 iter=0
869 if rtol>0:
870 r_dot_r = bilinearform(r, r)
871 if r_dot_r<0: raise NegativeNorm,"negative norm."
872 atol2=atol+rtol*math.sqrt(r_dot_r)
873 if verbose: print "GMRES: norm of right hand side = %e (absolute tolerance = %e)"%(math.sqrt(r_dot_r), atol2)
874 else:
875 atol2=atol
876 if verbose: print "GMRES: absolute tolerance = %e"%atol2
877 if atol2<=0:
878 raise ValueError,"Non-positive tolarance."
879
880 while True:
881 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached"%iter_max
882 if restarted:
883 r2 = r-Aprod(x-x2)
884 else:
885 r2=1*r
886 x2=x*1.
887 x,stopped=_GMRESm(r2, Aprod, x, bilinearform, atol2, iter_max=iter_max-iter, iter_restart=m, verbose=verbose,P_R=P_R)
888 iter+=iter_restart
889 if stopped: break
890 if verbose: print "GMRES: restart."
891 restarted=True
892 if verbose: print "GMRES: tolerance has been reached."
893 return x
894
895 def _GMRESm(r, Aprod, x, bilinearform, atol, iter_max=100, iter_restart=20, verbose=False, P_R=None):
896 iter=0
897
898 h=numarray.zeros((iter_restart+1,iter_restart),numarray.Float64)
899 c=numarray.zeros(iter_restart,numarray.Float64)
900 s=numarray.zeros(iter_restart,numarray.Float64)
901 g=numarray.zeros(iter_restart+1,numarray.Float64)
902 v=[]
903
904 r_dot_r = bilinearform(r, r)
905 if r_dot_r<0: raise NegativeNorm,"negative norm."
906 rho=math.sqrt(r_dot_r)
907
908 v.append(r/rho)
909 g[0]=rho
910
911 if verbose: print "GMRES: initial residual %e (absolute tolerance = %e)"%(rho,atol)
912 while not (rho<=atol or iter==iter_restart):
913
914 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
915
916 if P_R!=None:
917 p=Aprod(P_R(v[iter]))
918 else:
919 p=Aprod(v[iter])
920 v.append(p)
921
922 v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
923
924 # Modified Gram-Schmidt
925 for j in range(iter+1):
926 h[j,iter]=bilinearform(v[j],v[iter+1])
927 v[iter+1]-=h[j,iter]*v[j]
928
929 h[iter+1,iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
930 v_norm2=h[iter+1,iter]
931
932 # Reorthogonalize if needed
933 if v_norm1 + 0.001*v_norm2 == v_norm1: #Brown/Hindmarsh condition (default)
934 for j in range(iter+1):
935 hr=bilinearform(v[j],v[iter+1])
936 h[j,iter]=h[j,iter]+hr
937 v[iter+1] -= hr*v[j]
938
939 v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
940 h[iter+1,iter]=v_norm2
941
942 # watch out for happy breakdown
943 if not v_norm2 == 0:
944 v[iter+1]=v[iter+1]/h[iter+1,iter]
945
946 # Form and store the information for the new Givens rotation
947 if iter > 0: h[:iter+1,iter]=__givapp(c[:iter],s[:iter],h[:iter+1,iter])
948 mu=math.sqrt(h[iter,iter]*h[iter,iter]+h[iter+1,iter]*h[iter+1,iter])
949
950 if mu!=0 :
951 c[iter]=h[iter,iter]/mu
952 s[iter]=-h[iter+1,iter]/mu
953 h[iter,iter]=c[iter]*h[iter,iter]-s[iter]*h[iter+1,iter]
954 h[iter+1,iter]=0.0
955 gg=__givapp(c[iter],s[iter],[g[iter],g[iter+1]])
956 g[iter]=gg[0]
957 g[iter+1]=gg[1]
958 # Update the residual norm
959
960 rho=abs(g[iter+1])
961 if verbose: print "GMRES: iteration step %s: residual %e"%(iter,rho)
962 iter+=1
963
964 # At this point either iter > iter_max or rho < tol.
965 # It's time to compute x and leave.
966
967 if verbose: print "GMRES: iteration stopped after %s step."%iter
968 if iter > 0 :
969 y=numarray.zeros(iter,numarray.Float64)
970 y[iter-1] = g[iter-1] / h[iter-1,iter-1]
971 if iter > 1 :
972 i=iter-2
973 while i>=0 :
974 y[i] = ( g[i] - numarray.dot(h[i,i+1:iter], y[i+1:iter])) / h[i,i]
975 i=i-1
976 xhat=v[iter-1]*y[iter-1]
977 for i in range(iter-1):
978 xhat += v[i]*y[i]
979 else:
980 xhat=v[0] * 0
981 if P_R!=None:
982 x += P_R(xhat)
983 else:
984 x += xhat
985 if iter<iter_restart-1:
986 stopped=True
987 else:
988 stopped=False
989
990 return x,stopped
991
992 def MINRES(r, Aprod, x, Msolve, bilinearform, atol=0, rtol=1.e-8, iter_max=100):
993 """
994 Solver for
995
996 M{Ax=b}
997
998 with a symmetric and positive definite operator A (more details required!).
999 It uses the minimum residual method (MINRES) with preconditioner M
1000 providing an approximation of A.
1001
1002 The iteration is terminated if
1003
1004 M{|r| <= atol+rtol*|r0|}
1005
1006 where M{r0} is the initial residual and M{|.|} is the energy norm. In fact
1007
1008 M{|r| = sqrt( bilinearform(Msolve(r),r))}
1009
1010 For details on the preconditioned conjugate gradient method see the book:
1011
1012 I{Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
1013 T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
1014 C. Romine, and H. van der Vorst}.
1015
1016 @param r: initial residual M{r=b-Ax}. C{r} is altered.
1017 @type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1018 @param x: an initial guess for the solution
1019 @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1020 @param Aprod: returns the value Ax
1021 @type Aprod: function C{Aprod(x)} where C{x} is of the same object like
1022 argument C{x}. The returned object needs to be of the same
1023 type like argument C{r}.
1024 @param Msolve: solves Mx=r
1025 @type Msolve: function C{Msolve(r)} where C{r} is of the same type like
1026 argument C{r}. The returned object needs to be of the same
1027 type like argument C{x}.
1028 @param bilinearform: inner product C{<x,r>}
1029 @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same
1030 type like argument C{x} and C{r} is. The returned value
1031 is a C{float}.
1032 @param atol: absolute tolerance
1033 @type atol: non-negative C{float}
1034 @param rtol: relative tolerance
1035 @type rtol: non-negative C{float}
1036 @param iter_max: maximum number of iteration steps
1037 @type iter_max: C{int}
1038 @return: the solution approximation and the corresponding residual
1039 @rtype: C{tuple}
1040 @warning: C{r} and C{x} are altered.
1041 """
1042 #------------------------------------------------------------------
1043 # Set up y and v for the first Lanczos vector v1.
1044 # y = beta1 P' v1, where P = C**(-1).
1045 # v is really P' v1.
1046 #------------------------------------------------------------------
1047 r1 = r
1048 y = Msolve(r)
1049 beta1 = bilinearform(y,r)
1050
1051 if beta1< 0: raise NegativeNorm,"negative norm."
1052
1053 # If r = 0 exactly, stop with x
1054 if beta1==0: return x
1055
1056 if beta1> 0: beta1 = math.sqrt(beta1)
1057
1058 #------------------------------------------------------------------
1059 # Initialize quantities.
1060 # ------------------------------------------------------------------
1061 iter = 0
1062 Anorm = 0
1063 ynorm = 0
1064 oldb = 0
1065 beta = beta1
1066 dbar = 0
1067 epsln = 0
1068 phibar = beta1
1069 rhs1 = beta1
1070 rhs2 = 0
1071 rnorm = phibar
1072 tnorm2 = 0
1073 ynorm2 = 0
1074 cs = -1
1075 sn = 0
1076 w = r*0.
1077 w2 = r*0.
1078 r2 = r1
1079 eps = 0.0001
1080
1081 #---------------------------------------------------------------------
1082 # Main iteration loop.
1083 # --------------------------------------------------------------------
1084 while not rnorm<=atol+rtol*Anorm*ynorm: # checks ||r|| < (||A|| ||x||) * TOL
1085
1086 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
1087 iter = iter + 1
1088
1089 #-----------------------------------------------------------------
1090 # Obtain quantities for the next Lanczos vector vk+1, k = 1, 2,...
1091 # The general iteration is similar to the case k = 1 with v0 = 0:
1092 #
1093 # p1 = Operator * v1 - beta1 * v0,
1094 # alpha1 = v1'p1,
1095 # q2 = p2 - alpha1 * v1,
1096 # beta2^2 = q2'q2,
1097 # v2 = (1/beta2) q2.
1098 #
1099 # Again, y = betak P vk, where P = C**(-1).
1100 #-----------------------------------------------------------------
1101 s = 1/beta # Normalize previous vector (in y).
1102 v = s*y # v = vk if P = I
1103
1104 y = Aprod(v)
1105
1106 if iter >= 2:
1107 y = y - (beta/oldb)*r1
1108
1109 alfa = bilinearform(v,y) # alphak
1110 y += (- alfa/beta)*r2
1111 r1 = r2
1112 r2 = y
1113 y = Msolve(r2)
1114 oldb = beta # oldb = betak
1115 beta = bilinearform(y,r2) # beta = betak+1^2
1116 if beta < 0: raise NegativeNorm,"negative norm."
1117
1118 beta = math.sqrt( beta )
1119 tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta
1120
1121 if iter==1: # Initialize a few things.
1122 gmax = abs( alfa ) # alpha1
1123 gmin = gmax # alpha1
1124
1125 # Apply previous rotation Qk-1 to get
1126 # [deltak epslnk+1] = [cs sn][dbark 0 ]
1127 # [gbar k dbar k+1] [sn -cs][alfak betak+1].
1128
1129 oldeps = epsln
1130 delta = cs * dbar + sn * alfa # delta1 = 0 deltak
1131 gbar = sn * dbar - cs * alfa # gbar 1 = alfa1 gbar k
1132 epsln = sn * beta # epsln2 = 0 epslnk+1
1133 dbar = - cs * beta # dbar 2 = beta2 dbar k+1
1134
1135 # Compute the next plane rotation Qk
1136
1137 gamma = math.sqrt(gbar*gbar+beta*beta) # gammak
1138 gamma = max(gamma,eps)
1139 cs = gbar / gamma # ck
1140 sn = beta / gamma # sk
1141 phi = cs * phibar # phik
1142 phibar = sn * phibar # phibark+1
1143
1144 # Update x.
1145
1146 denom = 1/gamma
1147 w1 = w2
1148 w2 = w
1149 w = (v - oldeps*w1 - delta*w2) * denom
1150 x += phi*w
1151
1152 # Go round again.
1153
1154 gmax = max(gmax,gamma)
1155 gmin = min(gmin,gamma)
1156 z = rhs1 / gamma
1157 ynorm2 = z*z + ynorm2
1158 rhs1 = rhs2 - delta*z
1159 rhs2 = - epsln*z
1160
1161 # Estimate various norms and test for convergence.
1162
1163 Anorm = math.sqrt( tnorm2 )
1164 ynorm = math.sqrt( ynorm2 )
1165
1166 rnorm = phibar
1167
1168 return x
1169
1170 def TFQMR(r, Aprod, x, bilinearform, atol=0, rtol=1.e-8, iter_max=100):
1171 """
1172 Solver for
1173
1174 M{Ax=b}
1175
1176 with a general operator A (more details required!).
1177 It uses the Transpose-Free Quasi-Minimal Residual method (TFQMR).
1178
1179 The iteration is terminated if
1180
1181 M{|r| <= atol+rtol*|r0|}
1182
1183 where M{r0} is the initial residual and M{|.|} is the energy norm. In fact
1184
1185 M{|r| = sqrt( bilinearform(r,r))}
1186
1187 @param r: initial residual M{r=b-Ax}. C{r} is altered.
1188 @type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1189 @param x: an initial guess for the solution
1190 @type x: same like C{r}
1191 @param Aprod: returns the value Ax
1192 @type Aprod: function C{Aprod(x)} where C{x} is of the same object like
1193 argument C{x}. The returned object needs to be of the same type
1194 like argument C{r}.
1195 @param bilinearform: inner product C{<x,r>}
1196 @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same
1197 type like argument C{x} and C{r}. The returned value is
1198 a C{float}.
1199 @param atol: absolute tolerance
1200 @type atol: non-negative C{float}
1201 @param rtol: relative tolerance
1202 @type rtol: non-negative C{float}
1203 @param iter_max: maximum number of iteration steps
1204 @type iter_max: C{int}
1205 @rtype: C{tuple}
1206 @warning: C{r} and C{x} are altered.
1207 """
1208 u1=0
1209 u2=0
1210 y1=0
1211 y2=0
1212
1213 w = r
1214 y1 = r
1215 iter = 0
1216 d = 0
1217 v = Aprod(y1)
1218 u1 = v
1219
1220 theta = 0.0;
1221 eta = 0.0;
1222 rho=bilinearform(r,r)
1223 if rho < 0: raise NegativeNorm,"negative norm."
1224 tau = math.sqrt(rho)
1225 norm_r0=tau
1226 while tau>atol+rtol*norm_r0:
1227 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
1228
1229 sigma = bilinearform(r,v)
1230 if sigma == 0.0: raise IterationBreakDown,'TFQMR breakdown, sigma=0'
1231
1232 alpha = rho / sigma
1233
1234 for j in range(2):
1235 #
1236 # Compute y2 and u2 only if you have to
1237 #
1238 if ( j == 1 ):
1239 y2 = y1 - alpha * v
1240 u2 = Aprod(y2)
1241
1242 m = 2 * (iter+1) - 2 + (j+1)
1243 if j==0:
1244 w = w - alpha * u1
1245 d = y1 + ( theta * theta * eta / alpha ) * d
1246 if j==1:
1247 w = w - alpha * u2
1248 d = y2 + ( theta * theta * eta / alpha ) * d
1249
1250 theta = math.sqrt(bilinearform(w,w))/ tau
1251 c = 1.0 / math.sqrt ( 1.0 + theta * theta )
1252 tau = tau * theta * c
1253 eta = c * c * alpha
1254 x = x + eta * d
1255 #
1256 # Try to terminate the iteration at each pass through the loop
1257 #
1258 if rho == 0.0: raise IterationBreakDown,'TFQMR breakdown, rho=0'
1259
1260 rhon = bilinearform(r,w)
1261 beta = rhon / rho;
1262 rho = rhon;
1263 y1 = w + beta * y2;
1264 u1 = Aprod(y1)
1265 v = u1 + beta * ( u2 + beta * v )
1266
1267 iter += 1
1268
1269 return x
1270
1271
1272 #############################################
1273
1274 class ArithmeticTuple(object):
1275 """
1276 Tuple supporting inplace update x+=y and scaling x=a*y where C{x,y} is an
1277 ArithmeticTuple and C{a} is a float.
1278
1279 Example of usage::
1280
1281 from esys.escript import Data
1282 from numarray import array
1283 a=Data(...)
1284 b=array([1.,4.])
1285 x=ArithmeticTuple(a,b)
1286 y=5.*x
1287
1288 """
1289 def __init__(self,*args):
1290 """
1291 Initializes object with elements C{args}.
1292
1293 @param args: tuple of objects that support inplace add (x+=y) and
1294 scaling (x=a*y)
1295 """
1296 self.__items=list(args)
1297
1298 def __len__(self):
1299 """
1300 Returns the number of items.
1301
1302 @return: number of items
1303 @rtype: C{int}
1304 """
1305 return len(self.__items)
1306
1307 def __getitem__(self,index):
1308 """
1309 Returns item at specified position.
1310
1311 @param index: index of item to be returned
1312 @type index: C{int}
1313 @return: item with index C{index}
1314 """
1315 return self.__items.__getitem__(index)
1316
1317 def __mul__(self,other):
1318 """
1319 Scales by C{other} from the right.
1320
1321 @param other: scaling factor
1322 @type other: C{float}
1323 @return: itemwise self*other
1324 @rtype: L{ArithmeticTuple}
1325 """
1326 out=[]
1327 try:
1328 l=len(other)
1329 if l!=len(self):
1330 raise ValueError,"length of arguments don't match."
1331 for i in range(l): out.append(self[i]*other[i])
1332 except TypeError:
1333 for i in range(len(self)): out.append(self[i]*other)
1334 return ArithmeticTuple(*tuple(out))
1335
1336 def __rmul__(self,other):
1337 """
1338 Scales by C{other} from the left.
1339
1340 @param other: scaling factor
1341 @type other: C{float}
1342 @return: itemwise other*self
1343 @rtype: L{ArithmeticTuple}
1344 """
1345 out=[]
1346 try:
1347 l=len(other)
1348 if l!=len(self):
1349 raise ValueError,"length of arguments don't match."
1350 for i in range(l): out.append(other[i]*self[i])
1351 except TypeError:
1352 for i in range(len(self)): out.append(other*self[i])
1353 return ArithmeticTuple(*tuple(out))
1354
1355 def __div__(self,other):
1356 """
1357 Scales by (1/C{other}) from the right.
1358
1359 @param other: scaling factor
1360 @type other: C{float}
1361 @return: itemwise self/other
1362 @rtype: L{ArithmeticTuple}
1363 """
1364 return self*(1/other)
1365
1366 def __rdiv__(self,other):
1367 """
1368 Scales by (1/C{other}) from the left.
1369
1370 @param other: scaling factor
1371 @type other: C{float}
1372 @return: itemwise other/self
1373 @rtype: L{ArithmeticTuple}
1374 """
1375 out=[]
1376 try:
1377 l=len(other)
1378 if l!=len(self):
1379 raise ValueError,"length of arguments don't match."
1380 for i in range(l): out.append(other[i]/self[i])
1381 except TypeError:
1382 for i in range(len(self)): out.append(other/self[i])
1383 return ArithmeticTuple(*tuple(out))
1384
1385 def __iadd__(self,other):
1386 """
1387 Inplace addition of C{other} to self.
1388
1389 @param other: increment
1390 @type other: C{ArithmeticTuple}
1391 """
1392 if len(self) != len(other):
1393 raise ValueError,"tuple lengths must match."
1394 for i in range(len(self)):
1395 self.__items[i]+=other[i]
1396 return self
1397
1398 def __add__(self,other):
1399 """
1400 Adds C{other} to self.
1401
1402 @param other: increment
1403 @type other: C{ArithmeticTuple}
1404 """
1405 out=[]
1406 try:
1407 l=len(other)
1408 if l!=len(self):
1409 raise ValueError,"length of arguments don't match."
1410 for i in range(l): out.append(self[i]+other[i])
1411 except TypeError:
1412 for i in range(len(self)): out.append(self[i]+other)
1413 return ArithmeticTuple(*tuple(out))
1414
1415 def __sub__(self,other):
1416 """
1417 Subtracts C{other} from self.
1418
1419 @param other: decrement
1420 @type other: C{ArithmeticTuple}
1421 """
1422 out=[]
1423 try:
1424 l=len(other)
1425 if l!=len(self):
1426 raise ValueError,"length of arguments don't match."
1427 for i in range(l): out.append(self[i]-other[i])
1428 except TypeError:
1429 for i in range(len(self)): out.append(self[i]-other)
1430 return ArithmeticTuple(*tuple(out))
1431
1432 def __isub__(self,other):
1433 """
1434 Inplace subtraction of C{other} from self.
1435
1436 @param other: decrement
1437 @type other: C{ArithmeticTuple}
1438 """
1439 if len(self) != len(other):
1440 raise ValueError,"tuple length must match."
1441 for i in range(len(self)):
1442 self.__items[i]-=other[i]
1443 return self
1444
1445 def __neg__(self):
1446 """
1447 Negates values.
1448 """
1449 out=[]
1450 for i in range(len(self)):
1451 out.append(-self[i])
1452 return ArithmeticTuple(*tuple(out))
1453
1454
1455 class HomogeneousSaddlePointProblem(object):
1456 """
1457 This class provides a framework for solving linear homogeneous saddle
1458 point problems of the form::
1459
1460 M{Av+B^*p=f}
1461 M{Bv =0}
1462
1463 for the unknowns M{v} and M{p} and given operators M{A} and M{B} and
1464 given right hand side M{f}. M{B^*} is the adjoint operator of M{B}.
1465 """
1466 def __init__(self,**kwargs):
1467 self.setTolerance()
1468 self.setAbsoluteTolerance()
1469 self.setSubProblemTolerance()
1470
1471 #=============================================================
1472 def initialize(self):
1473 """
1474 Initializes the problem (overwrite).
1475 """
1476 pass
1477
1478 def inner_pBv(self,p,v):
1479 """
1480 Returns inner product of element p and Bv (overwrite).
1481
1482 @param p: a pressure increment
1483 @param v: a residual
1484 @return: inner product of element p and Bv
1485 @rtype: C{float}
1486 @note: used if PCG is applied.
1487 """
1488 raise NotImplementedError,"no inner product for p implemented."
1489
1490 def inner_p(self,p0,p1):
1491 """
1492 Returns inner product of p0 and p1 (overwrite).
1493
1494 @param p0: a pressure
1495 @param p1: a pressure
1496 @return: inner product of p0 and p1
1497 @rtype: C{float}
1498 """
1499 raise NotImplementedError,"no inner product for p implemented."
1500
1501 def norm_v(self,v):
1502 """
1503 Returns the norm of v (overwrite).
1504
1505 @param v: a velovity
1506 @return: norm of v
1507 @rtype: non-negative C{float}
1508 """
1509 raise NotImplementedError,"no norm of v implemented."
1510
1511
1512 def getV(self, p, v0):
1513 """
1514 return the value for v for a given p (overwrite)
1515
1516 @param p: a pressure
1517 @param v0: a initial guess for the value v to return.
1518 @return: v given as M{v= A^{-1} (f-B^*p)}
1519 """
1520 raise NotImplementedError,"no v calculation implemented."
1521
1522 def norm_Bv(self,v):
1523 """
1524 Returns Bv (overwrite).
1525
1526 @rtype: equal to the type of p
1527 @note: boundary conditions on p should be zero!
1528 """
1529 raise NotImplementedError, "no operator B implemented."
1530
1531 def solve_AinvBt(self,p):
1532 """
1533 Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}
1534 (overwrite).
1535
1536 @param p: a pressure increment
1537 @return: the solution of M{Av=B^*p}
1538 @note: boundary conditions on v should be zero!
1539 """
1540 raise NotImplementedError,"no operator A implemented."
1541
1542 def solve_precB(self,v):
1543 """
1544 Provides a preconditioner for M{BA^{-1}B^*} with accuracy
1545 L{self.getSubProblemTolerance()} (overwrite).
1546
1547 @rtype: equal to the type of p
1548 @note: boundary conditions on p should be zero!
1549 """
1550 raise NotImplementedError,"no preconditioner for Schur complement implemented."
1551 #=============================================================
1552 def __Aprod_PCG(self,p):
1553 return self.solve_AinvBt(p)
1554
1555 def __inner_PCG(self,p,v):
1556 return self.inner_pBv(p,v)
1557
1558 def __Msolve_PCG(self,v):
1559 return self.solve_precB(v)
1560 #=============================================================
1561 def __Aprod_GMRES(self,p):
1562 return self.solve_precB(self.solve_AinvBt(p))
1563 def __inner_GMRES(self,p0,p1):
1564 return self.inner_p(p0,p1)
1565 #=============================================================
1566 def norm_p(self,p):
1567 """
1568 calculates the norm of C{p}
1569
1570 @param p: a pressure
1571 @return: the norm of C{p} using the inner product for pressure
1572 @rtype: C{float}
1573 """
1574 f=self.inner_p(p,p)
1575 if f<0: raise ValueError,"negative pressure norm."
1576 return math.sqrt(f)
1577
1578
1579 def solve(self,v,p,max_iter=20, verbose=False, show_details=False, usePCG=True, iter_restart=20, max_correction_steps=10):
1580 """
1581 Solves the saddle point problem using initial guesses v and p.
1582
1583 @param v: initial guess for velocity
1584 @param p: initial guess for pressure
1585 @type v: L{Data}
1586 @type p: L{Data}
1587 @param usePCG: indicates the usage of the PCG rather than GMRES scheme.
1588 @param max_iter: maximum number of iteration steps per correction
1589 attempt
1590 @param verbose: if True, shows information on the progress of the
1591 saddlepoint problem solver.
1592 @param show_details: if True, shows details of the sub problem solver
1593 @param iter_restart: restart the iteration after C{iter_restart} steps
1594 (only used if useUzaw=False)
1595 @type usePCG: C{bool}
1596 @type max_iter: C{int}
1597 @type verbose: C{bool}
1598 @type show_details: C{bool}
1599 @type iter_restart: C{int}
1600 @rtype: C{tuple} of L{Data} objects
1601 """
1602 self.verbose=verbose
1603 self.show_details=show_details and self.verbose
1604 rtol=self.getTolerance()
1605 atol=self.getAbsoluteTolerance()
1606 correction_step=0
1607 converged=False
1608 while not converged:
1609 # calculate velocity for current pressure:
1610 v=self.getV(p,v)
1611 #
1612 norm_v=self.norm_v(v)
1613 norm_Bv=self.norm_Bv(v)
1614 ATOL=norm_v*rtol+atol
1615 if self.verbose: print "saddle point solver: norm v= %e, norm_Bv= %e, tolerance = %e."%(norm_v, norm_Bv,ATOL)
1616 if not ATOL>0: raise ValueError,"overall absolute tolerance needs to be positive."
1617 if norm_Bv <= ATOL:
1618 converged=True
1619 else:
1620 correction_step+=1
1621 if correction_step>max_correction_steps:
1622 raise CorrectionFailed,"Given up after %d correction steps."%correction_step
1623 dp=self.solve_precB(v)
1624 if usePCG:
1625 norm2=self.inner_pBv(dp,v)
1626 if norm2<0: raise ValueError,"negative PCG norm."
1627 norm2=math.sqrt(norm2)
1628 else:
1629 norm2=self.norm_p(dp)
1630 ATOL_ITER=ATOL/norm_Bv*norm2
1631 if self.verbose: print "saddle point solver: tolerance for solver: %e"%ATOL_ITER
1632 if usePCG:
1633 p,v0,a_norm=PCG(v,self.__Aprod_PCG,p,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL_ITER, rtol=0.,iter_max=max_iter, verbose=self.verbose)
1634 else:
1635 p=GMRES(dp,self.__Aprod_GMRES, p, self.__inner_GMRES,atol=ATOL_ITER, rtol=0.,iter_max=max_iter, iter_restart=iter_restart, verbose=self.verbose)
1636 if self.verbose: print "saddle point solver: tolerance reached."
1637 return v,p
1638
1639 #========================================================================
1640 def setTolerance(self,tolerance=1.e-4):
1641 """
1642 Sets the relative tolerance for (v,p).
1643
1644 @param tolerance: tolerance to be used
1645 @type tolerance: non-negative C{float}
1646 """
1647 if tolerance<0:
1648 raise ValueError,"tolerance must be positive."
1649 self.__rtol=tolerance
1650 self.setSubProblemTolerance()
1651
1652 def getTolerance(self):
1653 """
1654 Returns the relative tolerance.
1655
1656 @return: relative tolerance
1657 @rtype: C{float}
1658 """
1659 return self.__rtol
1660
1661 def setAbsoluteTolerance(self,tolerance=0.):
1662 """
1663 Sets the absolute tolerance.
1664
1665 @param tolerance: tolerance to be used
1666 @type tolerance: non-negative C{float}
1667 """
1668 if tolerance<0:
1669 raise ValueError,"tolerance must be non-negative."
1670 self.__atol=tolerance
1671
1672 def getAbsoluteTolerance(self):
1673 """
1674 Returns the absolute tolerance.
1675
1676 @return: absolute tolerance
1677 @rtype: C{float}
1678 """
1679 return self.__atol
1680
1681 def setSubProblemTolerance(self,rtol=None):
1682 """
1683 Sets the relative tolerance to solve the subproblem(s).
1684
1685 @param rtol: relative tolerence
1686 @type rtol: positive C{float}
1687 """
1688 if rtol == None:
1689 rtol=max(200.*util.EPSILON,self.getTolerance()**2)
1690 if rtol<=0:
1691 raise ValueError,"tolerance must be positive."
1692 self.__sub_tol=rtol
1693
1694 def getSubProblemTolerance(self):
1695 """
1696 Returns the subproblem reduction factor.
1697
1698 @return: subproblem reduction factor
1699 @rtype: C{float}
1700 """
1701 return self.__sub_tol
1702
1703 def MaskFromBoundaryTag(domain,*tags):
1704 """
1705 Creates a mask on the Solution(domain) function space where the value is
1706 one for samples that touch the boundary tagged by tags.
1707
1708 Usage: m=MaskFromBoundaryTag(domain, "left", "right")
1709
1710 @param domain: domain to be used
1711 @type domain: L{escript.Domain}
1712 @param tags: boundary tags
1713 @type tags: C{str}
1714 @return: a mask which marks samples that are touching the boundary tagged
1715 by any of the given tags
1716 @rtype: L{escript.Data} of rank 0
1717 """
1718 pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1719 d=escript.Scalar(0.,escript.FunctionOnBoundary(domain))
1720 for t in tags: d.setTaggedValue(t,1.)
1721 pde.setValue(y=d)
1722 return util.whereNonZero(pde.getRightHandSide())
1723
1724 #==============================================================================
1725 class SaddlePointProblem(object):
1726 """
1727 This implements a solver for a saddle point problem
1728
1729 M{f(u,p)=0}
1730 M{g(u)=0}
1731
1732 for u and p. The problem is solved with an inexact Uszawa scheme for p:
1733
1734 M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})}
1735 M{Q_g (p^{k+1}-p^{k}) = g(u^{k+1})}
1736
1737 where Q_f is an approximation of the Jacobian A_f of f with respect to u and
1738 Q_f is an approximation of A_g A_f^{-1} A_g with A_g is the Jacobian of g
1739 with respect to p. As a the construction of a 'proper' Q_g can be difficult,
1740 non-linear conjugate gradient method is applied to solve for p, so Q_g plays
1741 in fact the role of a preconditioner.
1742 """
1743 def __init__(self,verbose=False,*args):
1744 """
1745 Initializes the problem.
1746
1747 @param verbose: if True, some information is printed in the process
1748 @type verbose: C{bool}
1749 @note: this method may be overwritten by a particular saddle point
1750 problem
1751 """
1752 print "SaddlePointProblem should not be used anymore!"
1753 if not isinstance(verbose,bool):
1754 raise TypeError("verbose needs to be of type bool.")
1755 self.__verbose=verbose
1756 self.relaxation=1.
1757 DeprecationWarning("SaddlePointProblem should not be used anymore.")
1758
1759 def trace(self,text):
1760 """
1761 Prints C{text} if verbose has been set.
1762
1763 @param text: a text message
1764 @type text: C{str}
1765 """
1766 if self.__verbose: print "%s: %s"%(str(self),text)
1767
1768 def solve_f(self,u,p,tol=1.e-8):
1769 """
1770 Solves
1771
1772 A_f du = f(u,p)
1773
1774 with tolerance C{tol} and returns du. A_f is Jacobian of f with respect
1775 to u.
1776
1777 @param u: current approximation of u
1778 @type u: L{escript.Data}
1779 @param p: current approximation of p
1780 @type p: L{escript.Data}
1781 @param tol: tolerance expected for du
1782 @type tol: C{float}
1783 @return: increment du
1784 @rtype: L{escript.Data}
1785 @note: this method has to be overwritten by a particular saddle point
1786 problem
1787 """
1788 pass
1789
1790 def solve_g(self,u,tol=1.e-8):
1791 """
1792 Solves
1793
1794 Q_g dp = g(u)
1795
1796 where Q_g is a preconditioner for A_g A_f^{-1} A_g with A_g is the
1797 Jacobian of g with respect to p.
1798
1799 @param u: current approximation of u
1800 @type u: L{escript.Data}
1801 @param tol: tolerance expected for dp
1802 @type tol: C{float}
1803 @return: increment dp
1804 @rtype: L{escript.Data}
1805 @note: this method has to be overwritten by a particular saddle point
1806 problem
1807 """
1808 pass
1809
1810 def inner(self,p0,p1):
1811 """
1812 Inner product of p0 and p1 approximating p. Typically this returns
1813 C{integrate(p0*p1)}.
1814 @return: inner product of p0 and p1
1815 @rtype: C{float}
1816 """
1817 pass
1818
1819 subiter_max=3
1820 def solve(self,u0,p0,tolerance=1.e-6,tolerance_u=None,iter_max=100,accepted_reduction=0.995,relaxation=None):
1821 """
1822 Runs the solver.
1823
1824 @param u0: initial guess for C{u}
1825 @type u0: L{esys.escript.Data}
1826 @param p0: initial guess for C{p}
1827 @type p0: L{esys.escript.Data}
1828 @param tolerance: tolerance for relative error in C{u} and C{p}
1829 @type tolerance: positive C{float}
1830 @param tolerance_u: tolerance for relative error in C{u} if different
1831 from C{tolerance}
1832 @type tolerance_u: positive C{float}
1833 @param iter_max: maximum number of iteration steps
1834 @type iter_max: C{int}
1835 @param accepted_reduction: if the norm g cannot be reduced by
1836 C{accepted_reduction} backtracking to adapt
1837 the relaxation factor. If
1838 C{accepted_reduction=None} no backtracking
1839 is used.
1840 @type accepted_reduction: positive C{float} or C{None}
1841 @param relaxation: initial relaxation factor. If C{relaxation==None},
1842 the last relaxation factor is used.
1843 @type relaxation: C{float} or C{None}
1844 """
1845 tol=1.e-2
1846 if tolerance_u==None: tolerance_u=tolerance
1847 if not relaxation==None: self.relaxation=relaxation
1848 if accepted_reduction ==None:
1849 angle_limit=0.
1850 elif accepted_reduction>=1.:
1851 angle_limit=0.
1852 else:
1853 angle_limit=util.sqrt(1-accepted_reduction**2)
1854 self.iter=0
1855 u=u0
1856 p=p0
1857 #
1858 # initialize things:
1859 #
1860 converged=False
1861 #
1862 # start loop:
1863 #
1864 # initial search direction is g
1865 #
1866 while not converged :
1867 if self.iter>iter_max:
1868 raise ArithmeticError("no convergence after %s steps."%self.iter)
1869 f_new=self.solve_f(u,p,tol)
1870 norm_f_new = util.Lsup(f_new)
1871 u_new=u-f_new
1872 g_new=self.solve_g(u_new,tol)
1873 self.iter+=1
1874 norm_g_new = util.sqrt(self.inner(g_new,g_new))
1875 if norm_f_new==0. and norm_g_new==0.: return u, p
1876 if self.iter>1 and not accepted_reduction==None:
1877 #
1878 # did we manage to reduce the norm of G? I
1879 # if not we start a backtracking procedure
1880 #
1881 # print "new/old norm = ",norm_g_new, norm_g, norm_g_new/norm_g
1882 if norm_g_new > accepted_reduction * norm_g:
1883 sub_iter=0
1884 s=self.relaxation
1885 d=g
1886 g_last=g
1887 self.trace(" start substepping: f = %s, g = %s, relaxation = %s."%(norm_f_new, norm_g_new, s))
1888 while sub_iter < self.subiter_max and norm_g_new > accepted_reduction * norm_g:
1889 dg= g_new-g_last
1890 norm_dg=abs(util.sqrt(self.inner(dg,dg))/self.relaxation)
1891 rad=self.inner(g_new,dg)/self.relaxation
1892 # print " ",sub_iter,": rad, norm_dg:",abs(rad), norm_dg*norm_g_new * angle_limit
1893 # print " ",sub_iter,": rad, norm_dg:",rad, norm_dg, norm_g_new, norm_g
1894 if abs(rad) < norm_dg*norm_g_new * angle_limit:
1895 if sub_iter>0: self.trace(" no further improvements expected from backtracking.")
1896 break
1897 r=self.relaxation
1898 self.relaxation= - rad/norm_dg**2
1899 s+=self.relaxation
1900 #####
1901 # a=g_new+self.relaxation*dg/r
1902 # print "predicted new norm = ",util.sqrt(self.inner(a,a)),util.sqrt(self.inner(g_new,g_new)), self.relaxation
1903 #####
1904 g_last=g_new
1905 p+=self.relaxation*d
1906 f_new=self.solve_f(u,p,tol)
1907 u_new=u-f_new
1908 g_new=self.solve_g(u_new,tol)
1909 self.iter+=1
1910 norm_f_new = util.Lsup(f_new)
1911 norm_g_new = util.sqrt(self.inner(g_new,g_new))
1912 # print " ",sub_iter," new g norm",norm_g_new
1913 self.trace(" %s th sub-step: f = %s, g = %s, relaxation = %s."%(sub_iter, norm_f_new, norm_g_new, s))
1914 #
1915 # can we expect reduction of g?
1916 #
1917 # u_last=u_new
1918 sub_iter+=1
1919 self.relaxation=s
1920 #
1921 # check for convergence:
1922 #
1923 norm_u_new = util.Lsup(u_new)
1924 p_new=p+self.relaxation*g_new
1925 norm_p_new = util.sqrt(self.inner(p_new,p_new))
1926 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))
1927
1928 if self.iter>1:
1929 dg2=g_new-g
1930 df2=f_new-f
1931 norm_dg2=util.sqrt(self.inner(dg2,dg2))
1932 norm_df2=util.Lsup(df2)
1933 # print norm_g_new, norm_g, norm_dg, norm_p, tolerance
1934 tol_eq_g=tolerance*norm_dg2/(norm_g*abs(self.relaxation))*norm_p_new
1935 tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new
1936 if norm_g_new <= tol_eq_g and norm_f_new <= tol_eq_f:
1937 converged=True
1938 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
1939 self.trace("convergence after %s steps."%self.iter)
1940 return u,p
1941

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26