/[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 2455 - (show annotations)
Wed Jun 3 03:29:07 2009 UTC (10 years, 8 months ago) by jfenwick
File MIME type: text/x-python
File size: 64950 byte(s)
Merging changes from numpy branch.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26