/[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 2484 - (show annotations)
Mon Jun 22 04:22:19 2009 UTC (9 years, 6 months ago) by gross
File MIME type: text/x-python
File size: 56479 byte(s)
numarray removed from docu; Locator revised.
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__="https://launchpad.net/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 numpy
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.getSolverOptions().setSolverMethod(linearPDEs.SolverOptions.LUMPING)
149 self.__pde.setSymmetryOn()
150 self.__pde.setReducedOrderTo(reduce)
151 self.__pde.setValue(D = 1.)
152 return
153 def getSolverOptions(self):
154 """
155 Returns the solver options of the PDE solver.
156
157 @rtype: L{linearPDEs.SolverOptions}
158 """
159
160 def __call__(self, input_data):
161 """
162 Projects C{input_data} onto a continuous function.
163
164 @param input_data: the data to be projected
165 """
166 out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())
167 self.__pde.setValue(Y = escript.Data(), Y_reduced = escript.Data())
168 if input_data.getRank()==0:
169 self.__pde.setValue(Y = input_data)
170 out=self.__pde.getSolution()
171 elif input_data.getRank()==1:
172 for i0 in range(input_data.getShape()[0]):
173 self.__pde.setValue(Y = input_data[i0])
174 out[i0]=self.__pde.getSolution()
175 elif input_data.getRank()==2:
176 for i0 in range(input_data.getShape()[0]):
177 for i1 in range(input_data.getShape()[1]):
178 self.__pde.setValue(Y = input_data[i0,i1])
179 out[i0,i1]=self.__pde.getSolution()
180 elif input_data.getRank()==3:
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 self.__pde.setValue(Y = input_data[i0,i1,i2])
185 out[i0,i1,i2]=self.__pde.getSolution()
186 else:
187 for i0 in range(input_data.getShape()[0]):
188 for i1 in range(input_data.getShape()[1]):
189 for i2 in range(input_data.getShape()[2]):
190 for i3 in range(input_data.getShape()[3]):
191 self.__pde.setValue(Y = input_data[i0,i1,i2,i3])
192 out[i0,i1,i2,i3]=self.__pde.getSolution()
193 return out
194
195 class NoPDE:
196 """
197 Solves the following problem for u:
198
199 M{kronecker[i,j]*D[j]*u[j]=Y[i]}
200
201 with constraint
202
203 M{u[j]=r[j]} where M{q[j]>0}
204
205 where M{D}, M{Y}, M{r} and M{q} are given functions of rank 1.
206
207 In the case of scalars this takes the form
208
209 M{D*u=Y}
210
211 with constraint
212
213 M{u=r} where M{q>0}
214
215 where M{D}, M{Y}, M{r} and M{q} are given scalar functions.
216
217 The constraint overwrites any other condition.
218
219 @note: This class is similar to the L{linearPDEs.LinearPDE} class with
220 A=B=C=X=0 but has the intention that all input parameters are given
221 in L{Solution} or L{ReducedSolution}.
222 """
223 # The whole thing is a bit strange and I blame Rob Woodcock (CSIRO) for
224 # this.
225 def __init__(self,domain,D=None,Y=None,q=None,r=None):
226 """
227 Initializes the problem.
228
229 @param domain: domain of the PDE
230 @type domain: L{Domain}
231 @param D: coefficient of the solution
232 @type D: C{float}, C{int}, C{numpy.ndarray}, L{Data}
233 @param Y: right hand side
234 @type Y: C{float}, C{int}, C{numpy.ndarray}, L{Data}
235 @param q: location of constraints
236 @type q: C{float}, C{int}, C{numpy.ndarray}, L{Data}
237 @param r: value of solution at locations of constraints
238 @type r: C{float}, C{int}, C{numpy.ndarray}, L{Data}
239 """
240 self.__domain=domain
241 self.__D=D
242 self.__Y=Y
243 self.__q=q
244 self.__r=r
245 self.__u=None
246 self.__function_space=escript.Solution(self.__domain)
247
248 def setReducedOn(self):
249 """
250 Sets the L{FunctionSpace} of the solution to L{ReducedSolution}.
251 """
252 self.__function_space=escript.ReducedSolution(self.__domain)
253 self.__u=None
254
255 def setReducedOff(self):
256 """
257 Sets the L{FunctionSpace} of the solution to L{Solution}.
258 """
259 self.__function_space=escript.Solution(self.__domain)
260 self.__u=None
261
262 def setValue(self,D=None,Y=None,q=None,r=None):
263 """
264 Assigns values to the parameters.
265
266 @param D: coefficient of the solution
267 @type D: C{float}, C{int}, C{numpy.ndarray}, L{Data}
268 @param Y: right hand side
269 @type Y: C{float}, C{int}, C{numpy.ndarray}, L{Data}
270 @param q: location of constraints
271 @type q: C{float}, C{int}, C{numpy.ndarray}, L{Data}
272 @param r: value of solution at locations of constraints
273 @type r: C{float}, C{int}, C{numpy.ndarray}, L{Data}
274 """
275 if not D==None:
276 self.__D=D
277 self.__u=None
278 if not Y==None:
279 self.__Y=Y
280 self.__u=None
281 if not q==None:
282 self.__q=q
283 self.__u=None
284 if not r==None:
285 self.__r=r
286 self.__u=None
287
288 def getSolution(self):
289 """
290 Returns the solution.
291
292 @return: the solution of the problem
293 @rtype: L{Data} object in the L{FunctionSpace} L{Solution} or
294 L{ReducedSolution}
295 """
296 if self.__u==None:
297 if self.__D==None:
298 raise ValueError,"coefficient D is undefined"
299 D=escript.Data(self.__D,self.__function_space)
300 if D.getRank()>1:
301 raise ValueError,"coefficient D must have rank 0 or 1"
302 if self.__Y==None:
303 self.__u=escript.Data(0.,D.getShape(),self.__function_space)
304 else:
305 self.__u=util.quotient(self.__Y,D)
306 if not self.__q==None:
307 q=util.wherePositive(escript.Data(self.__q,self.__function_space))
308 self.__u*=(1.-q)
309 if not self.__r==None: self.__u+=q*self.__r
310 return self.__u
311
312 class Locator:
313 """
314 Locator provides access to the values of data objects at a given spatial
315 coordinate x.
316
317 In fact, a Locator object finds the sample in the set of samples of a
318 given function space or domain which is closest to the given point x.
319 """
320
321 def __init__(self,where,x=numpy.zeros((3,))):
322 """
323 Initializes a Locator to access values in Data objects on the Doamin
324 or FunctionSpace for the sample point which is closest to the given
325 point x.
326
327 @param where: function space
328 @type where: L{escript.FunctionSpace}
329 @param x: location(s) of the Locator
330 @type x: C{numpy.ndarray} or C{list} of C{numpy.ndarray}
331 """
332 if isinstance(where,escript.FunctionSpace):
333 self.__function_space=where
334 else:
335 self.__function_space=escript.ContinuousFunction(where)
336 iterative=False
337 if isinstance(x, list):
338 if len(x)==0:
339 raise "ValueError", "At least one point must be given."
340 try:
341 iter(x[0])
342 iterative=True
343 except TypeError:
344 iterative=False
345 if iterative:
346 self.__id=[]
347 for p in x:
348 self.__id.append(util.length(self.__function_space.getX()-p[:self.__function_space.getDim()]).minGlobalDataPoint())
349 else:
350 self.__id=util.length(self.__function_space.getX()-x[:self.__function_space.getDim()]).minGlobalDataPoint()
351
352 def __str__(self):
353 """
354 Returns the coordinates of the Locator as a string.
355 """
356 x=self.getX()
357 if instance(x,list):
358 out="["
359 first=True
360 for xx in x:
361 if not first:
362 out+=","
363 else:
364 first=False
365 out+=str(xx)
366 out+="]>"
367 else:
368 out=str(x)
369 return out
370
371 def getX(self):
372 """
373 Returns the exact coordinates of the Locator.
374 """
375 return self(self.getFunctionSpace().getX())
376
377 def getFunctionSpace(self):
378 """
379 Returns the function space of the Locator.
380 """
381 return self.__function_space
382
383 def getId(self,item=None):
384 """
385 Returns the identifier of the location.
386 """
387 if item == None:
388 return self.__id
389 else:
390 if isinstance(self.__id,list):
391 return self.__id[item]
392 else:
393 return self.__id
394
395
396 def __call__(self,data):
397 """
398 Returns the value of data at the Locator of a Data object.
399 """
400 return self.getValue(data)
401
402 def getValue(self,data):
403 """
404 Returns the value of C{data} at the Locator if C{data} is a L{Data}
405 object otherwise the object is returned.
406 """
407 if isinstance(data,escript.Data):
408 if data.getFunctionSpace()==self.getFunctionSpace():
409 dat=data
410 else:
411 dat=data.interpolate(self.getFunctionSpace())
412 id=self.getId()
413 r=data.getRank()
414 if isinstance(id,list):
415 out=[]
416 for i in id:
417 o=numpy.array(data.getTupleForGlobalDataPoint(*i))
418 if data.getRank()==0:
419 out.append(o[0])
420 else:
421 out.append(o)
422 return out
423 else:
424 out=numpy.array(data.getTupleForGlobalDataPoint(*id))
425 if data.getRank()==0:
426 return out[0]
427 else:
428 return out
429 else:
430 return data
431
432 class SolverSchemeException(Exception):
433 """
434 This is a generic exception thrown by solvers.
435 """
436 pass
437
438 class IndefinitePreconditioner(SolverSchemeException):
439 """
440 Exception thrown if the preconditioner is not positive definite.
441 """
442 pass
443
444 class MaxIterReached(SolverSchemeException):
445 """
446 Exception thrown if the maximum number of iteration steps is reached.
447 """
448 pass
449
450 class CorrectionFailed(SolverSchemeException):
451 """
452 Exception thrown if no convergence has been achieved in the solution
453 correction scheme.
454 """
455 pass
456
457 class IterationBreakDown(SolverSchemeException):
458 """
459 Exception thrown if the iteration scheme encountered an incurable breakdown.
460 """
461 pass
462
463 class NegativeNorm(SolverSchemeException):
464 """
465 Exception thrown if a norm calculation returns a negative norm.
466 """
467 pass
468
469 def PCG(r, Aprod, x, Msolve, bilinearform, atol=0, rtol=1.e-8, iter_max=100, initial_guess=True, verbose=False):
470 """
471 Solver for
472
473 M{Ax=b}
474
475 with a symmetric and positive definite operator A (more details required!).
476 It uses the conjugate gradient method with preconditioner M providing an
477 approximation of A.
478
479 The iteration is terminated if
480
481 M{|r| <= atol+rtol*|r0|}
482
483 where M{r0} is the initial residual and M{|.|} is the energy norm. In fact
484
485 M{|r| = sqrt( bilinearform(Msolve(r),r))}
486
487 For details on the preconditioned conjugate gradient method see the book:
488
489 I{Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
490 T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
491 C. Romine, and H. van der Vorst}.
492
493 @param r: initial residual M{r=b-Ax}. C{r} is altered.
494 @type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
495 @param x: an initial guess for the solution
496 @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
497 @param Aprod: returns the value Ax
498 @type Aprod: function C{Aprod(x)} where C{x} is of the same object like
499 argument C{x}. The returned object needs to be of the same type
500 like argument C{r}.
501 @param Msolve: solves Mx=r
502 @type Msolve: function C{Msolve(r)} where C{r} is of the same type like
503 argument C{r}. The returned object needs to be of the same
504 type like argument C{x}.
505 @param bilinearform: inner product C{<x,r>}
506 @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same
507 type like argument C{x} and C{r} is. The returned value
508 is a C{float}.
509 @param atol: absolute tolerance
510 @type atol: non-negative C{float}
511 @param rtol: relative tolerance
512 @type rtol: non-negative C{float}
513 @param iter_max: maximum number of iteration steps
514 @type iter_max: C{int}
515 @return: the solution approximation and the corresponding residual
516 @rtype: C{tuple}
517 @warning: C{r} and C{x} are altered.
518 """
519 iter=0
520 rhat=Msolve(r)
521 d = rhat
522 rhat_dot_r = bilinearform(rhat, r)
523 if rhat_dot_r<0: raise NegativeNorm,"negative norm."
524 norm_r0=math.sqrt(rhat_dot_r)
525 atol2=atol+rtol*norm_r0
526 if atol2<=0:
527 raise ValueError,"Non-positive tolarance."
528 atol2=max(atol2, 100. * util.EPSILON * norm_r0)
529
530 if verbose: print "PCG: initial residual norm = %e (absolute tolerance = %e)"%(norm_r0, atol2)
531
532
533 while not math.sqrt(rhat_dot_r) <= atol2:
534 iter+=1
535 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
536
537 q=Aprod(d)
538 alpha = rhat_dot_r / bilinearform(d, q)
539 x += alpha * d
540 if isinstance(q,ArithmeticTuple):
541 r += q * (-alpha) # Doing it the other way calls the float64.__mul__ not AT.__rmul__
542 else:
543 r += (-alpha) * q
544 rhat=Msolve(r)
545 rhat_dot_r_new = bilinearform(rhat, r)
546 beta = rhat_dot_r_new / rhat_dot_r
547 rhat+=beta * d
548 d=rhat
549
550 rhat_dot_r = rhat_dot_r_new
551 if rhat_dot_r<0: raise NegativeNorm,"negative norm."
552 if verbose: print "PCG: iteration step %s: residual norm = %e"%(iter, math.sqrt(rhat_dot_r))
553 if verbose: print "PCG: tolerance reached after %s steps."%iter
554 return x,r,math.sqrt(rhat_dot_r)
555
556 class Defect(object):
557 """
558 Defines a non-linear defect F(x) of a variable x.
559 """
560 def __init__(self):
561 """
562 Initializes defect.
563 """
564 self.setDerivativeIncrementLength()
565
566 def bilinearform(self, x0, x1):
567 """
568 Returns the inner product of x0 and x1
569
570 @param x0: value for x0
571 @param x1: value for x1
572 @return: the inner product of x0 and x1
573 @rtype: C{float}
574 """
575 return 0
576
577 def norm(self,x):
578 """
579 Returns the norm of argument C{x}.
580
581 @param x: a value
582 @return: norm of argument x
583 @rtype: C{float}
584 @note: by default C{sqrt(self.bilinearform(x,x)} is returned.
585 """
586 s=self.bilinearform(x,x)
587 if s<0: raise NegativeNorm,"negative norm."
588 return math.sqrt(s)
589
590 def eval(self,x):
591 """
592 Returns the value F of a given C{x}.
593
594 @param x: value for which the defect C{F} is evaluated
595 @return: value of the defect at C{x}
596 """
597 return 0
598
599 def __call__(self,x):
600 return self.eval(x)
601
602 def setDerivativeIncrementLength(self,inc=math.sqrt(util.EPSILON)):
603 """
604 Sets the relative length of the increment used to approximate the
605 derivative of the defect. The increment is inc*norm(x)/norm(v)*v in the
606 direction of v with x as a starting point.
607
608 @param inc: relative increment length
609 @type inc: positive C{float}
610 """
611 if inc<=0: raise ValueError,"positive increment required."
612 self.__inc=inc
613
614 def getDerivativeIncrementLength(self):
615 """
616 Returns the relative increment length used to approximate the
617 derivative of the defect.
618 @return: value of the defect at C{x}
619 @rtype: positive C{float}
620 """
621 return self.__inc
622
623 def derivative(self, F0, x0, v, v_is_normalised=True):
624 """
625 Returns the directional derivative at C{x0} in the direction of C{v}.
626
627 @param F0: value of this defect at x0
628 @param x0: value at which derivative is calculated
629 @param v: direction
630 @param v_is_normalised: True to indicate that C{v} is nomalized
631 (self.norm(v)=0)
632 @return: derivative of this defect at x0 in the direction of C{v}
633 @note: by default numerical evaluation (self.eval(x0+eps*v)-F0)/eps is
634 used but this method maybe overwritten to use exact evaluation.
635 """
636 normx=self.norm(x0)
637 if normx>0:
638 epsnew = self.getDerivativeIncrementLength() * normx
639 else:
640 epsnew = self.getDerivativeIncrementLength()
641 if not v_is_normalised:
642 normv=self.norm(v)
643 if normv<=0:
644 return F0*0
645 else:
646 epsnew /= normv
647 F1=self.eval(x0 + epsnew * v)
648 return (F1-F0)/epsnew
649
650 ######################################
651 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):
652 """
653 Solves a non-linear problem M{F(x)=0} for unknown M{x} using the stopping
654 criterion:
655
656 M{norm(F(x) <= atol + rtol * norm(F(x0)}
657
658 where M{x0} is the initial guess.
659
660 @param defect: object defining the function M{F}. C{defect.norm} defines the
661 M{norm} used in the stopping criterion.
662 @type defect: L{Defect}
663 @param x: initial guess for the solution, C{x} is altered.
664 @type x: any object type allowing basic operations such as
665 C{numpy.ndarray}, L{Data}
666 @param iter_max: maximum number of iteration steps
667 @type iter_max: positive C{int}
668 @param sub_iter_max: maximum number of inner iteration steps
669 @type sub_iter_max: positive C{int}
670 @param atol: absolute tolerance for the solution
671 @type atol: positive C{float}
672 @param rtol: relative tolerance for the solution
673 @type rtol: positive C{float}
674 @param gamma: tolerance safety factor for inner iteration
675 @type gamma: positive C{float}, less than 1
676 @param sub_tol_max: upper bound for inner tolerance
677 @type sub_tol_max: positive C{float}, less than 1
678 @return: an approximation of the solution with the desired accuracy
679 @rtype: same type as the initial guess
680 """
681 lmaxit=iter_max
682 if atol<0: raise ValueError,"atol needs to be non-negative."
683 if rtol<0: raise ValueError,"rtol needs to be non-negative."
684 if rtol+atol<=0: raise ValueError,"rtol or atol needs to be non-negative."
685 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
686 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
687
688 F=defect(x)
689 fnrm=defect.norm(F)
690 stop_tol=atol + rtol*fnrm
691 sub_tol=sub_tol_max
692 if verbose: print "NewtonGMRES: initial residual = %e."%fnrm
693 if verbose: print " tolerance = %e."%sub_tol
694 iter=1
695 #
696 # main iteration loop
697 #
698 while not fnrm<=stop_tol:
699 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
700 #
701 # adjust sub_tol_
702 #
703 if iter > 1:
704 rat=fnrm/fnrmo
705 sub_tol_old=sub_tol
706 sub_tol=gamma*rat**2
707 if gamma*sub_tol_old**2 > .1: sub_tol=max(sub_tol,gamma*sub_tol_old**2)
708 sub_tol=max(min(sub_tol,sub_tol_max), .5*stop_tol/fnrm)
709 #
710 # calculate newton increment xc
711 # if iter_max in __FDGMRES is reached MaxIterReached is thrown
712 # if iter_restart -1 is returned as sub_iter
713 # if atol is reached sub_iter returns the numer of steps performed to get there
714 #
715 #
716 if verbose: print " subiteration (GMRES) is called with relative tolerance %e."%sub_tol
717 try:
718 xc, sub_iter=__FDGMRES(F, defect, x, sub_tol*fnrm, iter_max=iter_max-iter, iter_restart=sub_iter_max)
719 except MaxIterReached:
720 raise MaxIterReached,"maximum number of %s steps reached."%iter_max
721 if sub_iter<0:
722 iter+=sub_iter_max
723 else:
724 iter+=sub_iter
725 # ====
726 x+=xc
727 F=defect(x)
728 iter+=1
729 fnrmo, fnrm=fnrm, defect.norm(F)
730 if verbose: print " step %s: residual %e."%(iter,fnrm)
731 if verbose: print "NewtonGMRES: completed after %s steps."%iter
732 return x
733
734 def __givapp(c,s,vin):
735 """
736 Applies a sequence of Givens rotations (c,s) recursively to the vector
737 C{vin}
738
739 @warning: C{vin} is altered.
740 """
741 vrot=vin
742 if isinstance(c,float):
743 vrot=[c*vrot[0]-s*vrot[1],s*vrot[0]+c*vrot[1]]
744 else:
745 for i in range(len(c)):
746 w1=c[i]*vrot[i]-s[i]*vrot[i+1]
747 w2=s[i]*vrot[i]+c[i]*vrot[i+1]
748 vrot[i]=w1
749 vrot[i+1]=w2
750 return vrot
751
752 def __FDGMRES(F0, defect, x0, atol, iter_max=100, iter_restart=20):
753 h=numpy.zeros((iter_restart,iter_restart),numpy.float64)
754 c=numpy.zeros(iter_restart,numpy.float64)
755 s=numpy.zeros(iter_restart,numpy.float64)
756 g=numpy.zeros(iter_restart,numpy.float64)
757 v=[]
758
759 rho=defect.norm(F0)
760 if rho<=0.: return x0*0
761
762 v.append(-F0/rho)
763 g[0]=rho
764 iter=0
765 while rho > atol and iter<iter_restart-1:
766 if iter >= iter_max:
767 raise MaxIterReached,"maximum number of %s steps reached."%iter_max
768
769 p=defect.derivative(F0,x0,v[iter], v_is_normalised=True)
770 v.append(p)
771
772 v_norm1=defect.norm(v[iter+1])
773
774 # Modified Gram-Schmidt
775 for j in range(iter+1):
776 h[j,iter]=defect.bilinearform(v[j],v[iter+1])
777 v[iter+1]-=h[j,iter]*v[j]
778
779 h[iter+1,iter]=defect.norm(v[iter+1])
780 v_norm2=h[iter+1,iter]
781
782 # Reorthogonalize if needed
783 if v_norm1 + 0.001*v_norm2 == v_norm1: #Brown/Hindmarsh condition (default)
784 for j in range(iter+1):
785 hr=defect.bilinearform(v[j],v[iter+1])
786 h[j,iter]=h[j,iter]+hr
787 v[iter+1] -= hr*v[j]
788
789 v_norm2=defect.norm(v[iter+1])
790 h[iter+1,iter]=v_norm2
791 # watch out for happy breakdown
792 if not v_norm2 == 0:
793 v[iter+1]=v[iter+1]/h[iter+1,iter]
794
795 # Form and store the information for the new Givens rotation
796 if iter > 0 :
797 hhat=numpy.zeros(iter+1,numpy.float64)
798 for i in range(iter+1) : hhat[i]=h[i,iter]
799 hhat=__givapp(c[0:iter],s[0:iter],hhat);
800 for i in range(iter+1) : h[i,iter]=hhat[i]
801
802 mu=math.sqrt(h[iter,iter]*h[iter,iter]+h[iter+1,iter]*h[iter+1,iter])
803
804 if mu!=0 :
805 c[iter]=h[iter,iter]/mu
806 s[iter]=-h[iter+1,iter]/mu
807 h[iter,iter]=c[iter]*h[iter,iter]-s[iter]*h[iter+1,iter]
808 h[iter+1,iter]=0.0
809 gg=__givapp(c[iter],s[iter],[g[iter],g[iter+1]])
810 g[iter]=gg[0]
811 g[iter+1]=gg[1]
812
813 # Update the residual norm
814 rho=abs(g[iter+1])
815 iter+=1
816
817 # At this point either iter > iter_max or rho < tol.
818 # It's time to compute x and leave.
819 if iter > 0 :
820 y=numpy.zeros(iter,numpy.float64)
821 y[iter-1] = g[iter-1] / h[iter-1,iter-1]
822 if iter > 1 :
823 i=iter-2
824 while i>=0 :
825 y[i] = ( g[i] - numpy.dot(h[i,i+1:iter], y[i+1:iter])) / h[i,i]
826 i=i-1
827 xhat=v[iter-1]*y[iter-1]
828 for i in range(iter-1):
829 xhat += v[i]*y[i]
830 else :
831 xhat=v[0] * 0
832
833 if iter<iter_restart-1:
834 stopped=iter
835 else:
836 stopped=-1
837
838 return xhat,stopped
839
840 def GMRES(r, Aprod, x, bilinearform, atol=0, rtol=1.e-8, iter_max=100, iter_restart=20, verbose=False,P_R=None):
841 """
842 Solver for
843
844 M{Ax=b}
845
846 with a general operator A (more details required!).
847 It uses the generalized minimum residual method (GMRES).
848
849 The iteration is terminated if
850
851 M{|r| <= atol+rtol*|r0|}
852
853 where M{r0} is the initial residual and M{|.|} is the energy norm. In fact
854
855 M{|r| = sqrt( bilinearform(r,r))}
856
857 @param r: initial residual M{r=b-Ax}. C{r} is altered.
858 @type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
859 @param x: an initial guess for the solution
860 @type x: same like C{r}
861 @param Aprod: returns the value Ax
862 @type Aprod: function C{Aprod(x)} where C{x} is of the same object like
863 argument C{x}. The returned object needs to be of the same
864 type like argument C{r}.
865 @param bilinearform: inner product C{<x,r>}
866 @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same
867 type like argument C{x} and C{r}. The returned value is
868 a C{float}.
869 @param atol: absolute tolerance
870 @type atol: non-negative C{float}
871 @param rtol: relative tolerance
872 @type rtol: non-negative C{float}
873 @param iter_max: maximum number of iteration steps
874 @type iter_max: C{int}
875 @param iter_restart: in order to save memory the orthogonalization process
876 is terminated after C{iter_restart} steps and the
877 iteration is restarted.
878 @type iter_restart: C{int}
879 @return: the solution approximation and the corresponding residual
880 @rtype: C{tuple}
881 @warning: C{r} and C{x} are altered.
882 """
883 m=iter_restart
884 restarted=False
885 iter=0
886 if rtol>0:
887 r_dot_r = bilinearform(r, r)
888 if r_dot_r<0: raise NegativeNorm,"negative norm."
889 atol2=atol+rtol*math.sqrt(r_dot_r)
890 if verbose: print "GMRES: norm of right hand side = %e (absolute tolerance = %e)"%(math.sqrt(r_dot_r), atol2)
891 else:
892 atol2=atol
893 if verbose: print "GMRES: absolute tolerance = %e"%atol2
894 if atol2<=0:
895 raise ValueError,"Non-positive tolarance."
896
897 while True:
898 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached"%iter_max
899 if restarted:
900 r2 = r-Aprod(x-x2)
901 else:
902 r2=1*r
903 x2=x*1.
904 x,stopped=_GMRESm(r2, Aprod, x, bilinearform, atol2, iter_max=iter_max-iter, iter_restart=m, verbose=verbose,P_R=P_R)
905 iter+=iter_restart
906 if stopped: break
907 if verbose: print "GMRES: restart."
908 restarted=True
909 if verbose: print "GMRES: tolerance has been reached."
910 return x
911
912 def _GMRESm(r, Aprod, x, bilinearform, atol, iter_max=100, iter_restart=20, verbose=False, P_R=None):
913 iter=0
914
915 h=numpy.zeros((iter_restart+1,iter_restart),numpy.float64)
916 c=numpy.zeros(iter_restart,numpy.float64)
917 s=numpy.zeros(iter_restart,numpy.float64)
918 g=numpy.zeros(iter_restart+1,numpy.float64)
919 v=[]
920
921 r_dot_r = bilinearform(r, r)
922 if r_dot_r<0: raise NegativeNorm,"negative norm."
923 rho=math.sqrt(r_dot_r)
924
925 v.append(r/rho)
926 g[0]=rho
927
928 if verbose: print "GMRES: initial residual %e (absolute tolerance = %e)"%(rho,atol)
929 while not (rho<=atol or iter==iter_restart):
930
931 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
932
933 if P_R!=None:
934 p=Aprod(P_R(v[iter]))
935 else:
936 p=Aprod(v[iter])
937 v.append(p)
938
939 v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
940
941 # Modified Gram-Schmidt
942 for j in range(iter+1):
943 h[j,iter]=bilinearform(v[j],v[iter+1])
944 v[iter+1]-=h[j,iter]*v[j]
945
946 h[iter+1,iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
947 v_norm2=h[iter+1,iter]
948
949 # Reorthogonalize if needed
950 if v_norm1 + 0.001*v_norm2 == v_norm1: #Brown/Hindmarsh condition (default)
951 for j in range(iter+1):
952 hr=bilinearform(v[j],v[iter+1])
953 h[j,iter]=h[j,iter]+hr
954 v[iter+1] -= hr*v[j]
955
956 v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
957 h[iter+1,iter]=v_norm2
958
959 # watch out for happy breakdown
960 if not v_norm2 == 0:
961 v[iter+1]=v[iter+1]/h[iter+1,iter]
962
963 # Form and store the information for the new Givens rotation
964 if iter > 0: h[:iter+1,iter]=__givapp(c[:iter],s[:iter],h[:iter+1,iter])
965 mu=math.sqrt(h[iter,iter]*h[iter,iter]+h[iter+1,iter]*h[iter+1,iter])
966
967 if mu!=0 :
968 c[iter]=h[iter,iter]/mu
969 s[iter]=-h[iter+1,iter]/mu
970 h[iter,iter]=c[iter]*h[iter,iter]-s[iter]*h[iter+1,iter]
971 h[iter+1,iter]=0.0
972 gg=__givapp(c[iter],s[iter],[g[iter],g[iter+1]])
973 g[iter]=gg[0]
974 g[iter+1]=gg[1]
975 # Update the residual norm
976
977 rho=abs(g[iter+1])
978 if verbose: print "GMRES: iteration step %s: residual %e"%(iter,rho)
979 iter+=1
980
981 # At this point either iter > iter_max or rho < tol.
982 # It's time to compute x and leave.
983
984 if verbose: print "GMRES: iteration stopped after %s step."%iter
985 if iter > 0 :
986 y=numpy.zeros(iter,numpy.float64)
987 y[iter-1] = g[iter-1] / h[iter-1,iter-1]
988 if iter > 1 :
989 i=iter-2
990 while i>=0 :
991 y[i] = ( g[i] - numpy.dot(h[i,i+1:iter], y[i+1:iter])) / h[i,i]
992 i=i-1
993 xhat=v[iter-1]*y[iter-1]
994 for i in range(iter-1):
995 xhat += v[i]*y[i]
996 else:
997 xhat=v[0] * 0
998 if P_R!=None:
999 x += P_R(xhat)
1000 else:
1001 x += xhat
1002 if iter<iter_restart-1:
1003 stopped=True
1004 else:
1005 stopped=False
1006
1007 return x,stopped
1008
1009 def MINRES(r, Aprod, x, Msolve, bilinearform, atol=0, rtol=1.e-8, iter_max=100):
1010 """
1011 Solver for
1012
1013 M{Ax=b}
1014
1015 with a symmetric and positive definite operator A (more details required!).
1016 It uses the minimum residual method (MINRES) with preconditioner M
1017 providing an approximation of A.
1018
1019 The iteration is terminated if
1020
1021 M{|r| <= atol+rtol*|r0|}
1022
1023 where M{r0} is the initial residual and M{|.|} is the energy norm. In fact
1024
1025 M{|r| = sqrt( bilinearform(Msolve(r),r))}
1026
1027 For details on the preconditioned conjugate gradient method see the book:
1028
1029 I{Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
1030 T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
1031 C. Romine, and H. van der Vorst}.
1032
1033 @param r: initial residual M{r=b-Ax}. C{r} is altered.
1034 @type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1035 @param x: an initial guess for the solution
1036 @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1037 @param Aprod: returns the value Ax
1038 @type Aprod: function C{Aprod(x)} where C{x} is of the same object like
1039 argument C{x}. The returned object needs to be of the same
1040 type like argument C{r}.
1041 @param Msolve: solves Mx=r
1042 @type Msolve: function C{Msolve(r)} where C{r} is of the same type like
1043 argument C{r}. The returned object needs to be of the same
1044 type like argument C{x}.
1045 @param bilinearform: inner product C{<x,r>}
1046 @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same
1047 type like argument C{x} and C{r} is. The returned value
1048 is a C{float}.
1049 @param atol: absolute tolerance
1050 @type atol: non-negative C{float}
1051 @param rtol: relative tolerance
1052 @type rtol: non-negative C{float}
1053 @param iter_max: maximum number of iteration steps
1054 @type iter_max: C{int}
1055 @return: the solution approximation and the corresponding residual
1056 @rtype: C{tuple}
1057 @warning: C{r} and C{x} are altered.
1058 """
1059 #------------------------------------------------------------------
1060 # Set up y and v for the first Lanczos vector v1.
1061 # y = beta1 P' v1, where P = C**(-1).
1062 # v is really P' v1.
1063 #------------------------------------------------------------------
1064 r1 = r
1065 y = Msolve(r)
1066 beta1 = bilinearform(y,r)
1067
1068 if beta1< 0: raise NegativeNorm,"negative norm."
1069
1070 # If r = 0 exactly, stop with x
1071 if beta1==0: return x
1072
1073 if beta1> 0: beta1 = math.sqrt(beta1)
1074
1075 #------------------------------------------------------------------
1076 # Initialize quantities.
1077 # ------------------------------------------------------------------
1078 iter = 0
1079 Anorm = 0
1080 ynorm = 0
1081 oldb = 0
1082 beta = beta1
1083 dbar = 0
1084 epsln = 0
1085 phibar = beta1
1086 rhs1 = beta1
1087 rhs2 = 0
1088 rnorm = phibar
1089 tnorm2 = 0
1090 ynorm2 = 0
1091 cs = -1
1092 sn = 0
1093 w = r*0.
1094 w2 = r*0.
1095 r2 = r1
1096 eps = 0.0001
1097
1098 #---------------------------------------------------------------------
1099 # Main iteration loop.
1100 # --------------------------------------------------------------------
1101 while not rnorm<=atol+rtol*Anorm*ynorm: # checks ||r|| < (||A|| ||x||) * TOL
1102
1103 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
1104 iter = iter + 1
1105
1106 #-----------------------------------------------------------------
1107 # Obtain quantities for the next Lanczos vector vk+1, k = 1, 2,...
1108 # The general iteration is similar to the case k = 1 with v0 = 0:
1109 #
1110 # p1 = Operator * v1 - beta1 * v0,
1111 # alpha1 = v1'p1,
1112 # q2 = p2 - alpha1 * v1,
1113 # beta2^2 = q2'q2,
1114 # v2 = (1/beta2) q2.
1115 #
1116 # Again, y = betak P vk, where P = C**(-1).
1117 #-----------------------------------------------------------------
1118 s = 1/beta # Normalize previous vector (in y).
1119 v = s*y # v = vk if P = I
1120
1121 y = Aprod(v)
1122
1123 if iter >= 2:
1124 y = y - (beta/oldb)*r1
1125
1126 alfa = bilinearform(v,y) # alphak
1127 y += (- alfa/beta)*r2
1128 r1 = r2
1129 r2 = y
1130 y = Msolve(r2)
1131 oldb = beta # oldb = betak
1132 beta = bilinearform(y,r2) # beta = betak+1^2
1133 if beta < 0: raise NegativeNorm,"negative norm."
1134
1135 beta = math.sqrt( beta )
1136 tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta
1137
1138 if iter==1: # Initialize a few things.
1139 gmax = abs( alfa ) # alpha1
1140 gmin = gmax # alpha1
1141
1142 # Apply previous rotation Qk-1 to get
1143 # [deltak epslnk+1] = [cs sn][dbark 0 ]
1144 # [gbar k dbar k+1] [sn -cs][alfak betak+1].
1145
1146 oldeps = epsln
1147 delta = cs * dbar + sn * alfa # delta1 = 0 deltak
1148 gbar = sn * dbar - cs * alfa # gbar 1 = alfa1 gbar k
1149 epsln = sn * beta # epsln2 = 0 epslnk+1
1150 dbar = - cs * beta # dbar 2 = beta2 dbar k+1
1151
1152 # Compute the next plane rotation Qk
1153
1154 gamma = math.sqrt(gbar*gbar+beta*beta) # gammak
1155 gamma = max(gamma,eps)
1156 cs = gbar / gamma # ck
1157 sn = beta / gamma # sk
1158 phi = cs * phibar # phik
1159 phibar = sn * phibar # phibark+1
1160
1161 # Update x.
1162
1163 denom = 1/gamma
1164 w1 = w2
1165 w2 = w
1166 w = (v - oldeps*w1 - delta*w2) * denom
1167 x += phi*w
1168
1169 # Go round again.
1170
1171 gmax = max(gmax,gamma)
1172 gmin = min(gmin,gamma)
1173 z = rhs1 / gamma
1174 ynorm2 = z*z + ynorm2
1175 rhs1 = rhs2 - delta*z
1176 rhs2 = - epsln*z
1177
1178 # Estimate various norms and test for convergence.
1179
1180 Anorm = math.sqrt( tnorm2 )
1181 ynorm = math.sqrt( ynorm2 )
1182
1183 rnorm = phibar
1184
1185 return x
1186
1187 def TFQMR(r, Aprod, x, bilinearform, atol=0, rtol=1.e-8, iter_max=100):
1188 """
1189 Solver for
1190
1191 M{Ax=b}
1192
1193 with a general operator A (more details required!).
1194 It uses the Transpose-Free Quasi-Minimal Residual method (TFQMR).
1195
1196 The iteration is terminated if
1197
1198 M{|r| <= atol+rtol*|r0|}
1199
1200 where M{r0} is the initial residual and M{|.|} is the energy norm. In fact
1201
1202 M{|r| = sqrt( bilinearform(r,r))}
1203
1204 @param r: initial residual M{r=b-Ax}. C{r} is altered.
1205 @type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1206 @param x: an initial guess for the solution
1207 @type x: same like C{r}
1208 @param Aprod: returns the value Ax
1209 @type Aprod: function C{Aprod(x)} where C{x} is of the same object like
1210 argument C{x}. The returned object needs to be of the same type
1211 like argument C{r}.
1212 @param bilinearform: inner product C{<x,r>}
1213 @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same
1214 type like argument C{x} and C{r}. The returned value is
1215 a C{float}.
1216 @param atol: absolute tolerance
1217 @type atol: non-negative C{float}
1218 @param rtol: relative tolerance
1219 @type rtol: non-negative C{float}
1220 @param iter_max: maximum number of iteration steps
1221 @type iter_max: C{int}
1222 @rtype: C{tuple}
1223 @warning: C{r} and C{x} are altered.
1224 """
1225 u1=0
1226 u2=0
1227 y1=0
1228 y2=0
1229
1230 w = r
1231 y1 = r
1232 iter = 0
1233 d = 0
1234 v = Aprod(y1)
1235 u1 = v
1236
1237 theta = 0.0;
1238 eta = 0.0;
1239 rho=bilinearform(r,r)
1240 if rho < 0: raise NegativeNorm,"negative norm."
1241 tau = math.sqrt(rho)
1242 norm_r0=tau
1243 while tau>atol+rtol*norm_r0:
1244 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
1245
1246 sigma = bilinearform(r,v)
1247 if sigma == 0.0: raise IterationBreakDown,'TFQMR breakdown, sigma=0'
1248
1249 alpha = rho / sigma
1250
1251 for j in range(2):
1252 #
1253 # Compute y2 and u2 only if you have to
1254 #
1255 if ( j == 1 ):
1256 y2 = y1 - alpha * v
1257 u2 = Aprod(y2)
1258
1259 m = 2 * (iter+1) - 2 + (j+1)
1260 if j==0:
1261 w = w - alpha * u1
1262 d = y1 + ( theta * theta * eta / alpha ) * d
1263 if j==1:
1264 w = w - alpha * u2
1265 d = y2 + ( theta * theta * eta / alpha ) * d
1266
1267 theta = math.sqrt(bilinearform(w,w))/ tau
1268 c = 1.0 / math.sqrt ( 1.0 + theta * theta )
1269 tau = tau * theta * c
1270 eta = c * c * alpha
1271 x = x + eta * d
1272 #
1273 # Try to terminate the iteration at each pass through the loop
1274 #
1275 if rho == 0.0: raise IterationBreakDown,'TFQMR breakdown, rho=0'
1276
1277 rhon = bilinearform(r,w)
1278 beta = rhon / rho;
1279 rho = rhon;
1280 y1 = w + beta * y2;
1281 u1 = Aprod(y1)
1282 v = u1 + beta * ( u2 + beta * v )
1283
1284 iter += 1
1285
1286 return x
1287
1288
1289 #############################################
1290
1291 class ArithmeticTuple(object):
1292 """
1293 Tuple supporting inplace update x+=y and scaling x=a*y where C{x,y} is an
1294 ArithmeticTuple and C{a} is a float.
1295
1296 Example of usage::
1297
1298 from esys.escript import Data
1299 from numpy import array
1300 a=Data(...)
1301 b=array([1.,4.])
1302 x=ArithmeticTuple(a,b)
1303 y=5.*x
1304
1305 """
1306 def __init__(self,*args):
1307 """
1308 Initializes object with elements C{args}.
1309
1310 @param args: tuple of objects that support inplace add (x+=y) and
1311 scaling (x=a*y)
1312 """
1313 self.__items=list(args)
1314
1315 def __len__(self):
1316 """
1317 Returns the number of items.
1318
1319 @return: number of items
1320 @rtype: C{int}
1321 """
1322 return len(self.__items)
1323
1324 def __getitem__(self,index):
1325 """
1326 Returns item at specified position.
1327
1328 @param index: index of item to be returned
1329 @type index: C{int}
1330 @return: item with index C{index}
1331 """
1332 return self.__items.__getitem__(index)
1333
1334 def __mul__(self,other):
1335 """
1336 Scales by C{other} from the right.
1337
1338 @param other: scaling factor
1339 @type other: C{float}
1340 @return: itemwise self*other
1341 @rtype: L{ArithmeticTuple}
1342 """
1343 out=[]
1344 try:
1345 l=len(other)
1346 if l!=len(self):
1347 raise ValueError,"length of arguments don't match."
1348 for i in range(l): out.append(self[i]*other[i])
1349 except TypeError:
1350 for i in range(len(self)): out.append(self[i]*other)
1351 return ArithmeticTuple(*tuple(out))
1352
1353 def __rmul__(self,other):
1354 """
1355 Scales by C{other} from the left.
1356
1357 @param other: scaling factor
1358 @type other: C{float}
1359 @return: itemwise other*self
1360 @rtype: L{ArithmeticTuple}
1361 """
1362 out=[]
1363 try:
1364 l=len(other)
1365 if l!=len(self):
1366 raise ValueError,"length of arguments don't match."
1367 for i in range(l): out.append(other[i]*self[i])
1368 except TypeError:
1369 for i in range(len(self)): out.append(other*self[i])
1370 return ArithmeticTuple(*tuple(out))
1371
1372 def __div__(self,other):
1373 """
1374 Scales by (1/C{other}) from the right.
1375
1376 @param other: scaling factor
1377 @type other: C{float}
1378 @return: itemwise self/other
1379 @rtype: L{ArithmeticTuple}
1380 """
1381 return self*(1/other)
1382
1383 def __rdiv__(self,other):
1384 """
1385 Scales by (1/C{other}) from the left.
1386
1387 @param other: scaling factor
1388 @type other: C{float}
1389 @return: itemwise other/self
1390 @rtype: L{ArithmeticTuple}
1391 """
1392 out=[]
1393 try:
1394 l=len(other)
1395 if l!=len(self):
1396 raise ValueError,"length of arguments don't match."
1397 for i in range(l): out.append(other[i]/self[i])
1398 except TypeError:
1399 for i in range(len(self)): out.append(other/self[i])
1400 return ArithmeticTuple(*tuple(out))
1401
1402 def __iadd__(self,other):
1403 """
1404 Inplace addition of C{other} to self.
1405
1406 @param other: increment
1407 @type other: C{ArithmeticTuple}
1408 """
1409 if len(self) != len(other):
1410 raise ValueError,"tuple lengths must match."
1411 for i in range(len(self)):
1412 self.__items[i]+=other[i]
1413 return self
1414
1415 def __add__(self,other):
1416 """
1417 Adds C{other} to self.
1418
1419 @param other: increment
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 __sub__(self,other):
1433 """
1434 Subtracts C{other} from self.
1435
1436 @param other: decrement
1437 @type other: C{ArithmeticTuple}
1438 """
1439 out=[]
1440 try:
1441 l=len(other)
1442 if l!=len(self):
1443 raise ValueError,"length of arguments don't match."
1444 for i in range(l): out.append(self[i]-other[i])
1445 except TypeError:
1446 for i in range(len(self)): out.append(self[i]-other)
1447 return ArithmeticTuple(*tuple(out))
1448
1449 def __isub__(self,other):
1450 """
1451 Inplace subtraction of C{other} from self.
1452
1453 @param other: decrement
1454 @type other: C{ArithmeticTuple}
1455 """
1456 if len(self) != len(other):
1457 raise ValueError,"tuple length must match."
1458 for i in range(len(self)):
1459 self.__items[i]-=other[i]
1460 return self
1461
1462 def __neg__(self):
1463 """
1464 Negates values.
1465 """
1466 out=[]
1467 for i in range(len(self)):
1468 out.append(-self[i])
1469 return ArithmeticTuple(*tuple(out))
1470
1471
1472 class HomogeneousSaddlePointProblem(object):
1473 """
1474 This class provides a framework for solving linear homogeneous saddle
1475 point problems of the form::
1476
1477 M{Av+B^*p=f}
1478 M{Bv =0}
1479
1480 for the unknowns M{v} and M{p} and given operators M{A} and M{B} and
1481 given right hand side M{f}. M{B^*} is the adjoint operator of M{B}.
1482 """
1483 def __init__(self, adaptSubTolerance=True, **kwargs):
1484 """
1485 initializes the saddle point problem
1486
1487 @param adaptSubTolerance: If True the tolerance for subproblem is set automatically.
1488 @type adaptSubTolerance: C{bool}
1489 """
1490 self.setTolerance()
1491 self.setAbsoluteTolerance()
1492 self.__adaptSubTolerance=adaptSubTolerance
1493 #=============================================================
1494 def initialize(self):
1495 """
1496 Initializes the problem (overwrite).
1497 """
1498 pass
1499
1500 def inner_pBv(self,p,Bv):
1501 """
1502 Returns inner product of element p and Bv (overwrite).
1503
1504 @param p: a pressure increment
1505 @param v: a residual
1506 @return: inner product of element p and Bv
1507 @rtype: C{float}
1508 @note: used if PCG is applied.
1509 """
1510 raise NotImplementedError,"no inner product for p and Bv implemented."
1511
1512 def inner_p(self,p0,p1):
1513 """
1514 Returns inner product of p0 and p1 (overwrite).
1515
1516 @param p0: a pressure
1517 @param p1: a pressure
1518 @return: inner product of p0 and p1
1519 @rtype: C{float}
1520 """
1521 raise NotImplementedError,"no inner product for p implemented."
1522
1523 def norm_v(self,v):
1524 """
1525 Returns the norm of v (overwrite).
1526
1527 @param v: a velovity
1528 @return: norm of v
1529 @rtype: non-negative C{float}
1530 """
1531 raise NotImplementedError,"no norm of v implemented."
1532 def getV(self, p, v0):
1533 """
1534 return the value for v for a given p (overwrite)
1535
1536 @param p: a pressure
1537 @param v0: a initial guess for the value v to return.
1538 @return: v given as M{v= A^{-1} (f-B^*p)}
1539 """
1540 raise NotImplementedError,"no v calculation implemented."
1541
1542
1543 def Bv(self,v):
1544 """
1545 Returns Bv (overwrite).
1546
1547 @rtype: equal to the type of p
1548 @note: boundary conditions on p should be zero!
1549 """
1550 raise NotImplementedError, "no operator B implemented."
1551
1552 def norm_Bv(self,Bv):
1553 """
1554 Returns the norm of Bv (overwrite).
1555
1556 @rtype: equal to the type of p
1557 @note: boundary conditions on p should be zero!
1558 """
1559 raise NotImplementedError, "no norm of Bv implemented."
1560
1561 def solve_AinvBt(self,p):
1562 """
1563 Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}
1564 (overwrite).
1565
1566 @param p: a pressure increment
1567 @return: the solution of M{Av=B^*p}
1568 @note: boundary conditions on v should be zero!
1569 """
1570 raise NotImplementedError,"no operator A implemented."
1571
1572 def solve_prec(self,Bv):
1573 """
1574 Provides a preconditioner for M{BA^{-1}B^*} applied to Bv with accuracy
1575 L{self.getSubProblemTolerance()} (overwrite).
1576
1577 @rtype: equal to the type of p
1578 @note: boundary conditions on p should be zero!
1579 """
1580 raise NotImplementedError,"no preconditioner for Schur complement implemented."
1581 def setSubProblemTolerance(self):
1582 """
1583 Updates the tolerance for subproblems
1584 @note: method is typically the method is overwritten.
1585 """
1586 pass
1587 #=============================================================
1588 def __Aprod_PCG(self,p):
1589 dv=self.solve_AinvBt(p)
1590 return ArithmeticTuple(dv,self.Bv(dv))
1591
1592 def __inner_PCG(self,p,r):
1593 return self.inner_pBv(p,r[1])
1594
1595 def __Msolve_PCG(self,r):
1596 return self.solve_prec(r[1])
1597 #=============================================================
1598 def __Aprod_GMRES(self,p):
1599 return self.solve_prec(self.Bv(self.solve_AinvBt(p)))
1600 def __inner_GMRES(self,p0,p1):
1601 return self.inner_p(p0,p1)
1602
1603 #=============================================================
1604 def norm_p(self,p):
1605 """
1606 calculates the norm of C{p}
1607
1608 @param p: a pressure
1609 @return: the norm of C{p} using the inner product for pressure
1610 @rtype: C{float}
1611 """
1612 f=self.inner_p(p,p)
1613 if f<0: raise ValueError,"negative pressure norm."
1614 return math.sqrt(f)
1615 def adaptSubTolerance(self):
1616 """
1617 Returns True if tolerance adaption for subproblem is choosen.
1618 """
1619 self.__adaptSubTolerance
1620
1621 def solve(self,v,p,max_iter=20, verbose=False, usePCG=True, iter_restart=20, max_correction_steps=10):
1622 """
1623 Solves the saddle point problem using initial guesses v and p.
1624
1625 @param v: initial guess for velocity
1626 @param p: initial guess for pressure
1627 @type v: L{Data}
1628 @type p: L{Data}
1629 @param usePCG: indicates the usage of the PCG rather than GMRES scheme.
1630 @param max_iter: maximum number of iteration steps per correction
1631 attempt
1632 @param verbose: if True, shows information on the progress of the
1633 saddlepoint problem solver.
1634 @param iter_restart: restart the iteration after C{iter_restart} steps
1635 (only used if useUzaw=False)
1636 @type usePCG: C{bool}
1637 @type max_iter: C{int}
1638 @type verbose: C{bool}
1639 @type iter_restart: C{int}
1640 @rtype: C{tuple} of L{Data} objects
1641 """
1642 self.verbose=verbose
1643 rtol=self.getTolerance()
1644 atol=self.getAbsoluteTolerance()
1645 if self.adaptSubTolerance(): self.setSubProblemTolerance()
1646 correction_step=0
1647 converged=False
1648 while not converged:
1649 # calculate velocity for current pressure:
1650 v=self.getV(p,v)
1651 Bv=self.Bv(v)
1652 norm_v=self.norm_v(v)
1653 norm_Bv=self.norm_Bv(Bv)
1654 ATOL=norm_v*rtol+atol
1655 if self.verbose: print "HomogeneousSaddlePointProblem: norm v= %e, norm_Bv= %e, tolerance = %e."%(norm_v, norm_Bv,ATOL)
1656 if not ATOL>0: raise ValueError,"overall absolute tolerance needs to be positive."
1657 if norm_Bv <= ATOL:
1658 converged=True
1659 else:
1660 correction_step+=1
1661 if correction_step>max_correction_steps:
1662 raise CorrectionFailed,"Given up after %d correction steps."%correction_step
1663 dp=self.solve_prec(Bv)
1664 if usePCG:
1665 norm2=self.inner_pBv(dp,Bv)
1666 if norm2<0: raise ValueError,"negative PCG norm."
1667 norm2=math.sqrt(norm2)
1668 else:
1669 norm2=self.norm_p(dp)
1670 ATOL_ITER=ATOL/norm_Bv*norm2*0.5
1671 if self.verbose: print "HomogeneousSaddlePointProblem: tolerance for solver: %e"%ATOL_ITER
1672 if usePCG:
1673 p,v0,a_norm=PCG(ArithmeticTuple(v,Bv),self.__Aprod_PCG,p,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL_ITER, rtol=0.,iter_max=max_iter, verbose=self.verbose)
1674 else:
1675 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)
1676 if self.verbose: print "HomogeneousSaddlePointProblem: tolerance reached."
1677 return v,p
1678
1679 #========================================================================
1680 def setTolerance(self,tolerance=1.e-4):
1681 """
1682 Sets the relative tolerance for (v,p).
1683
1684 @param tolerance: tolerance to be used
1685 @type tolerance: non-negative C{float}
1686 """
1687 if tolerance<0:
1688 raise ValueError,"tolerance must be positive."
1689 self.__rtol=tolerance
1690
1691 def getTolerance(self):
1692 """
1693 Returns the relative tolerance.
1694
1695 @return: relative tolerance
1696 @rtype: C{float}
1697 """
1698 return self.__rtol
1699
1700 def setAbsoluteTolerance(self,tolerance=0.):
1701 """
1702 Sets the absolute tolerance.
1703
1704 @param tolerance: tolerance to be used
1705 @type tolerance: non-negative C{float}
1706 """
1707 if tolerance<0:
1708 raise ValueError,"tolerance must be non-negative."
1709 self.__atol=tolerance
1710
1711 def getAbsoluteTolerance(self):
1712 """
1713 Returns the absolute tolerance.
1714
1715 @return: absolute tolerance
1716 @rtype: C{float}
1717 """
1718 return self.__atol
1719
1720 def getSubProblemTolerance(self):
1721 """
1722 Sets the relative tolerance to solve the subproblem(s).
1723
1724 @param rtol: relative tolerence
1725 @type rtol: positive C{float}
1726 """
1727 return max(200.*util.EPSILON,self.getTolerance()**2)
1728
1729 def MaskFromBoundaryTag(domain,*tags):
1730 """
1731 Creates a mask on the Solution(domain) function space where the value is
1732 one for samples that touch the boundary tagged by tags.
1733
1734 Usage: m=MaskFromBoundaryTag(domain, "left", "right")
1735
1736 @param domain: domain to be used
1737 @type domain: L{escript.Domain}
1738 @param tags: boundary tags
1739 @type tags: C{str}
1740 @return: a mask which marks samples that are touching the boundary tagged
1741 by any of the given tags
1742 @rtype: L{escript.Data} of rank 0
1743 """
1744 pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1745 d=escript.Scalar(0.,escript.FunctionOnBoundary(domain))
1746 for t in tags: d.setTaggedValue(t,1.)
1747 pde.setValue(y=d)
1748 return util.whereNonZero(pde.getRightHandSide())
1749
1750

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26