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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26