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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26