/[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 2100 - (show annotations)
Wed Nov 26 08:13:00 2008 UTC (10 years, 4 months ago) by gross
File MIME type: text/x-python
File size: 62936 byte(s)
This commit cleans up the incompressible solver and adds a DarcyFlux solver in model module. 
Some documentation for both classes has been added.
The convection code is only linear at the moment.



1
2 ########################################################
3 #
4 # Copyright (c) 2003-2008 by University of Queensland
5 # Earth Systems Science Computational Center (ESSCC)
6 # http://www.uq.edu.au/esscc
7 #
8 # Primary Business: Queensland, Australia
9 # Licensed under the Open Software License version 3.0
10 # http://www.opensource.org/licenses/osl-3.0.php
11 #
12 ########################################################
13
14 __copyright__="""Copyright (c) 2003-2008 by University of Queensland
15 Earth Systems Science Computational Center (ESSCC)
16 http://www.uq.edu.au/esscc
17 Primary Business: Queensland, Australia"""
18 __license__="""Licensed under the Open Software License version 3.0
19 http://www.opensource.org/licenses/osl-3.0.php"""
20 __url__="http://www.uq.edu.au/esscc/escript-finley"
21
22 """
23 Provides some tools related to PDEs.
24
25 Currently includes:
26 - Projector - to project a discontinuous
27 - Locator - to trace values in data objects at a certain location
28 - TimeIntegrationManager - to handel extraplotion in time
29 - SaddlePointProblem - solver for Saddle point problems using the inexact uszawa scheme
30
31 @var __author__: name of author
32 @var __copyright__: copyrights
33 @var __license__: licence agreement
34 @var __url__: url entry point on documentation
35 @var __version__: version
36 @var __date__: date of the version
37 """
38
39 __author__="Lutz Gross, l.gross@uq.edu.au"
40
41
42 import escript
43 import linearPDEs
44 import numarray
45 import util
46 import math
47
48 ##### Added by Artak
49 # from Numeric import zeros,Int,Float64
50 ###################################
51
52
53 class TimeIntegrationManager:
54 """
55 a simple mechanism to manage time dependend values.
56
57 typical usage is::
58
59 dt=0.1 # time increment
60 tm=TimeIntegrationManager(inital_value,p=1)
61 while t<1.
62 v_guess=tm.extrapolate(dt) # extrapolate to t+dt
63 v=...
64 tm.checkin(dt,v)
65 t+=dt
66
67 @note: currently only p=1 is supported.
68 """
69 def __init__(self,*inital_values,**kwargs):
70 """
71 sets up the value manager where inital_value is the initial value and p is order used for extrapolation
72 """
73 if kwargs.has_key("p"):
74 self.__p=kwargs["p"]
75 else:
76 self.__p=1
77 if kwargs.has_key("time"):
78 self.__t=kwargs["time"]
79 else:
80 self.__t=0.
81 self.__v_mem=[inital_values]
82 self.__order=0
83 self.__dt_mem=[]
84 self.__num_val=len(inital_values)
85
86 def getTime(self):
87 return self.__t
88 def getValue(self):
89 out=self.__v_mem[0]
90 if len(out)==1:
91 return out[0]
92 else:
93 return out
94
95 def checkin(self,dt,*values):
96 """
97 adds new values to the manager. the p+1 last value get lost
98 """
99 o=min(self.__order+1,self.__p)
100 self.__order=min(self.__order+1,self.__p)
101 v_mem_new=[values]
102 dt_mem_new=[dt]
103 for i in range(o-1):
104 v_mem_new.append(self.__v_mem[i])
105 dt_mem_new.append(self.__dt_mem[i])
106 v_mem_new.append(self.__v_mem[o-1])
107 self.__order=o
108 self.__v_mem=v_mem_new
109 self.__dt_mem=dt_mem_new
110 self.__t+=dt
111
112 def extrapolate(self,dt):
113 """
114 extrapolates to dt forward in time.
115 """
116 if self.__order==0:
117 out=self.__v_mem[0]
118 else:
119 out=[]
120 for i in range(self.__num_val):
121 out.append((1.+dt/self.__dt_mem[0])*self.__v_mem[0][i]-dt/self.__dt_mem[0]*self.__v_mem[1][i])
122
123 if len(out)==0:
124 return None
125 elif len(out)==1:
126 return out[0]
127 else:
128 return out
129
130
131 class Projector:
132 """
133 The Projector is a factory which projects a discontiuous function onto a
134 continuous function on the a given domain.
135 """
136 def __init__(self, domain, reduce = True, fast=True):
137 """
138 Create a continuous function space projector for a domain.
139
140 @param domain: Domain of the projection.
141 @param reduce: Flag to reduce projection order (default is True)
142 @param fast: Flag to use a fast method based on matrix lumping (default is true)
143 """
144 self.__pde = linearPDEs.LinearPDE(domain)
145 if fast:
146 self.__pde.setSolverMethod(linearPDEs.LinearPDE.LUMPING)
147 self.__pde.setSymmetryOn()
148 self.__pde.setReducedOrderTo(reduce)
149 self.__pde.setValue(D = 1.)
150 return
151
152 def __call__(self, input_data):
153 """
154 Projects input_data onto a continuous function
155
156 @param input_data: The input_data to be projected.
157 """
158 out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())
159 self.__pde.setValue(Y = escript.Data(), Y_reduced = escript.Data())
160 if input_data.getRank()==0:
161 self.__pde.setValue(Y = input_data)
162 out=self.__pde.getSolution()
163 elif input_data.getRank()==1:
164 for i0 in range(input_data.getShape()[0]):
165 self.__pde.setValue(Y = input_data[i0])
166 out[i0]=self.__pde.getSolution()
167 elif input_data.getRank()==2:
168 for i0 in range(input_data.getShape()[0]):
169 for i1 in range(input_data.getShape()[1]):
170 self.__pde.setValue(Y = input_data[i0,i1])
171 out[i0,i1]=self.__pde.getSolution()
172 elif input_data.getRank()==3:
173 for i0 in range(input_data.getShape()[0]):
174 for i1 in range(input_data.getShape()[1]):
175 for i2 in range(input_data.getShape()[2]):
176 self.__pde.setValue(Y = input_data[i0,i1,i2])
177 out[i0,i1,i2]=self.__pde.getSolution()
178 else:
179 for i0 in range(input_data.getShape()[0]):
180 for i1 in range(input_data.getShape()[1]):
181 for i2 in range(input_data.getShape()[2]):
182 for i3 in range(input_data.getShape()[3]):
183 self.__pde.setValue(Y = input_data[i0,i1,i2,i3])
184 out[i0,i1,i2,i3]=self.__pde.getSolution()
185 return out
186
187 class NoPDE:
188 """
189 solves the following problem for u:
190
191 M{kronecker[i,j]*D[j]*u[j]=Y[i]}
192
193 with constraint
194
195 M{u[j]=r[j]} where M{q[j]>0}
196
197 where D, Y, r and q are given functions of rank 1.
198
199 In the case of scalars this takes the form
200
201 M{D*u=Y}
202
203 with constraint
204
205 M{u=r} where M{q>0}
206
207 where D, Y, r and q are given scalar functions.
208
209 The constraint is overwriting any other condition.
210
211 @note: This class is similar to the L{linearPDEs.LinearPDE} class with A=B=C=X=0 but has the intention
212 that all input parameter are given in L{Solution} or L{ReducedSolution}. The whole
213 thing is a bit strange and I blame Robert.Woodcock@csiro.au for this.
214 """
215 def __init__(self,domain,D=None,Y=None,q=None,r=None):
216 """
217 initialize the problem
218
219 @param domain: domain of the PDE.
220 @type domain: L{Domain}
221 @param D: coefficient of the solution.
222 @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}
223 @param Y: right hand side
224 @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}
225 @param q: location of constraints
226 @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}
227 @param r: value of solution at locations of constraints
228 @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}
229 """
230 self.__domain=domain
231 self.__D=D
232 self.__Y=Y
233 self.__q=q
234 self.__r=r
235 self.__u=None
236 self.__function_space=escript.Solution(self.__domain)
237 def setReducedOn(self):
238 """
239 sets the L{FunctionSpace} of the solution to L{ReducedSolution}
240 """
241 self.__function_space=escript.ReducedSolution(self.__domain)
242 self.__u=None
243
244 def setReducedOff(self):
245 """
246 sets the L{FunctionSpace} of the solution to L{Solution}
247 """
248 self.__function_space=escript.Solution(self.__domain)
249 self.__u=None
250
251 def setValue(self,D=None,Y=None,q=None,r=None):
252 """
253 assigns values to the parameters.
254
255 @param D: coefficient of the solution.
256 @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}
257 @param Y: right hand side
258 @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}
259 @param q: location of constraints
260 @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}
261 @param r: value of solution at locations of constraints
262 @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}
263 """
264 if not D==None:
265 self.__D=D
266 self.__u=None
267 if not Y==None:
268 self.__Y=Y
269 self.__u=None
270 if not q==None:
271 self.__q=q
272 self.__u=None
273 if not r==None:
274 self.__r=r
275 self.__u=None
276
277 def getSolution(self):
278 """
279 returns the solution
280
281 @return: the solution of the problem
282 @rtype: L{Data} object in the L{FunctionSpace} L{Solution} or L{ReducedSolution}.
283 """
284 if self.__u==None:
285 if self.__D==None:
286 raise ValueError,"coefficient D is undefined"
287 D=escript.Data(self.__D,self.__function_space)
288 if D.getRank()>1:
289 raise ValueError,"coefficient D must have rank 0 or 1"
290 if self.__Y==None:
291 self.__u=escript.Data(0.,D.getShape(),self.__function_space)
292 else:
293 self.__u=util.quotient(self.__Y,D)
294 if not self.__q==None:
295 q=util.wherePositive(escript.Data(self.__q,self.__function_space))
296 self.__u*=(1.-q)
297 if not self.__r==None: self.__u+=q*self.__r
298 return self.__u
299
300 class Locator:
301 """
302 Locator provides access to the values of data objects at a given
303 spatial coordinate x.
304
305 In fact, a Locator object finds the sample in the set of samples of a
306 given function space or domain where which is closest to the given
307 point x.
308 """
309
310 def __init__(self,where,x=numarray.zeros((3,))):
311 """
312 Initializes a Locator to access values in Data objects on the Doamin
313 or FunctionSpace where for the sample point which
314 closest to the given point x.
315
316 @param where: function space
317 @type where: L{escript.FunctionSpace}
318 @param x: coefficient of the solution.
319 @type x: L{numarray.NumArray} or C{list} of L{numarray.NumArray}
320 """
321 if isinstance(where,escript.FunctionSpace):
322 self.__function_space=where
323 else:
324 self.__function_space=escript.ContinuousFunction(where)
325 if isinstance(x, list):
326 self.__id=[]
327 for p in x:
328 self.__id.append(util.length(self.__function_space.getX()-p[:self.__function_space.getDim()]).minGlobalDataPoint())
329 else:
330 self.__id=util.length(self.__function_space.getX()-x[:self.__function_space.getDim()]).minGlobalDataPoint()
331
332 def __str__(self):
333 """
334 Returns the coordinates of the Locator as a string.
335 """
336 x=self.getX()
337 if instance(x,list):
338 out="["
339 first=True
340 for xx in x:
341 if not first:
342 out+=","
343 else:
344 first=False
345 out+=str(xx)
346 out+="]>"
347 else:
348 out=str(x)
349 return out
350
351 def getX(self):
352 """
353 Returns the exact coordinates of the Locator.
354 """
355 return self(self.getFunctionSpace().getX())
356
357 def getFunctionSpace(self):
358 """
359 Returns the function space of the Locator.
360 """
361 return self.__function_space
362
363 def getId(self,item=None):
364 """
365 Returns the identifier of the location.
366 """
367 if item == None:
368 return self.__id
369 else:
370 if isinstance(self.__id,list):
371 return self.__id[item]
372 else:
373 return self.__id
374
375
376 def __call__(self,data):
377 """
378 Returns the value of data at the Locator of a Data object otherwise
379 the object is returned.
380 """
381 return self.getValue(data)
382
383 def getValue(self,data):
384 """
385 Returns the value of data at the Locator if data is a Data object
386 otherwise the object is returned.
387 """
388 if isinstance(data,escript.Data):
389 if data.getFunctionSpace()==self.getFunctionSpace():
390 dat=data
391 else:
392 dat=data.interpolate(self.getFunctionSpace())
393 id=self.getId()
394 r=data.getRank()
395 if isinstance(id,list):
396 out=[]
397 for i in id:
398 o=data.getValueOfGlobalDataPoint(*i)
399 if data.getRank()==0:
400 out.append(o[0])
401 else:
402 out.append(o)
403 return out
404 else:
405 out=data.getValueOfGlobalDataPoint(*id)
406 if data.getRank()==0:
407 return out[0]
408 else:
409 return out
410 else:
411 return data
412
413 class SolverSchemeException(Exception):
414 """
415 exceptions thrown by solvers
416 """
417 pass
418
419 class IndefinitePreconditioner(SolverSchemeException):
420 """
421 the preconditioner is not positive definite.
422 """
423 pass
424 class MaxIterReached(SolverSchemeException):
425 """
426 maxium number of iteration steps is reached.
427 """
428 pass
429 class CorrectionFailed(SolverSchemeException):
430 """
431 no convergence has been achieved in the solution correction scheme.
432 """
433 pass
434 class IterationBreakDown(SolverSchemeException):
435 """
436 iteration scheme econouters an incurable breakdown.
437 """
438 pass
439 class NegativeNorm(SolverSchemeException):
440 """
441 a norm calculation returns a negative norm.
442 """
443 pass
444
445 def PCG(b, Aprod, Msolve, bilinearform, atol=0, rtol=1.e-8, x=None, iter_max=100, initial_guess=True, verbose=False):
446 """
447 Solver for
448
449 M{Ax=b}
450
451 with a symmetric and positive definite operator A (more details required!).
452 It uses the conjugate gradient method with preconditioner M providing an approximation of A.
453
454 The iteration is terminated if the C{stoppingcriterium} function return C{True}.
455
456 For details on the preconditioned conjugate gradient method see the book:
457
458 Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
459 T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
460 C. Romine, and H. van der Vorst.
461
462 @param b: the right hand side of the liner system. C{b} is altered.
463 @type b: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
464 @param Aprod: returns the value Ax
465 @type Aprod: function C{Aprod(x)} where C{x} is of the same object like argument C{x}.
466 The returned object needs to be of the same type like argument C{b}.
467 @param Msolve: solves Mx=r
468 @type Msolve: function C{Msolve(r)} where C{r} is of the same type like argument C{b}.
469 The returned object needs to be of the same type like argument C{x}.
470 @param bilinearform: inner product C{<x,r>}
471 @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same type like argument C{x} and C{r} is.
472 The returned value is a C{float}.
473 @param stoppingcriterium: function which returns True if a stopping criterium is meet.
474 C{stoppingcriterium} has the arguments C{norm_r}, C{r} and C{x} giving the current
475 norm of the residual (=C{sqrt(bilinearform(Msolve(r),r)}), the current residual and
476 the current solution approximation. C{stoppingcriterium} is called in each iteration step.
477 @type stoppingcriterium: function that returns C{True} or C{False}
478 @param x: an initial guess for the solution. If no C{x} is given 0*b is used.
479 @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
480 @param iter_max: maximum number of iteration steps.
481 @type iter_max: C{int}
482 @return: the solution approximation and the corresponding residual
483 @rtype: C{tuple}
484 @warning: C{b} and C{x} are altered.
485 """
486 iter=0
487 if x==None:
488 x=0*b
489 r=b*1.
490 initial_guess=False
491 if initial_guess:
492 r=b+(-1)*Aprod(x)
493 x=x*1.
494 else:
495 x=x*0
496 r=b*1
497 rhat=Msolve(r)
498 d = rhat
499 rhat_dot_r = bilinearform(rhat, r)
500 if rhat_dot_r<0: raise NegativeNorm,"negative norm."
501 norm_b=math.sqrt(rhat_dot_r)
502 if verbose: print "PCG: initial residual norm = %e (absolute tolerance = %e)"%(norm_b, atol+rtol*norm_b)
503
504 while not math.sqrt(rhat_dot_r) <= atol+rtol*norm_b:
505 iter+=1
506 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
507
508 q=Aprod(d)
509 alpha = rhat_dot_r / bilinearform(d, q)
510 x += alpha * d
511 r += (-alpha) * q
512
513 rhat=Msolve(r)
514 rhat_dot_r_new = bilinearform(rhat, r)
515 beta = rhat_dot_r_new / rhat_dot_r
516 rhat+=beta * d
517 d=rhat
518
519 rhat_dot_r = rhat_dot_r_new
520 if rhat_dot_r<0: raise NegativeNorm,"negative norm."
521 if verbose: print "PCG: iteration step %s: residual norm = %e"%(iter, math.sqrt(rhat_dot_r))
522 if verbose: print "PCG: tolerance reached after %s steps."%iter
523 return x,r
524
525 class Defect(object):
526 """
527 defines a non-linear defect F(x) of a variable x
528 """
529 def __init__(self):
530 """
531 initialize defect
532 """
533 self.setDerivativeIncrementLength()
534
535 def bilinearform(self, x0, x1):
536 """
537 returns the inner product of x0 and x1
538 @param x0: a value for x
539 @param x1: a value for x
540 @return: the inner product of x0 and x1
541 @rtype: C{float}
542 """
543 return 0
544
545 def norm(self,x):
546 """
547 the norm of argument C{x}
548
549 @param x: a value for x
550 @return: norm of argument x
551 @rtype: C{float}
552 @note: by default C{sqrt(self.bilinearform(x,x)} is retrurned.
553 """
554 s=self.bilinearform(x,x)
555 if s<0: raise NegativeNorm,"negative norm."
556 return math.sqrt(s)
557
558
559 def eval(self,x):
560 """
561 returns the value F of a given x
562
563 @param x: value for which the defect C{F} is evalulated.
564 @return: value of the defect at C{x}
565 """
566 return 0
567
568 def __call__(self,x):
569 return self.eval(x)
570
571 def setDerivativeIncrementLength(self,inc=math.sqrt(util.EPSILON)):
572 """
573 sets the relative length of the increment used to approximate the derivative of the defect
574 the increment is inc*norm(x)/norm(v)*v in the direction of v with x as a staring point.
575
576 @param inc: relative increment length
577 @type inc: positive C{float}
578 """
579 if inc<=0: raise ValueError,"positive increment required."
580 self.__inc=inc
581
582 def getDerivativeIncrementLength(self):
583 """
584 returns the relative increment length used to approximate the derivative of the defect
585 @return: value of the defect at C{x}
586 @rtype: positive C{float}
587 """
588 return self.__inc
589
590 def derivative(self, F0, x0, v, v_is_normalised=True):
591 """
592 returns the directional derivative at x0 in the direction of v
593
594 @param F0: value of this defect at x0
595 @param x0: value at which derivative is calculated.
596 @param v: direction
597 @param v_is_normalised: is true to indicate that C{v} is nomalized (self.norm(v)=0)
598 @return: derivative of this defect at x0 in the direction of C{v}
599 @note: by default numerical evaluation (self.eval(x0+eps*v)-F0)/eps is used but this method
600 maybe oepsnew verwritten to use exact evalution.
601 """
602 normx=self.norm(x0)
603 if normx>0:
604 epsnew = self.getDerivativeIncrementLength() * normx
605 else:
606 epsnew = self.getDerivativeIncrementLength()
607 if not v_is_normalised:
608 normv=self.norm(v)
609 if normv<=0:
610 return F0*0
611 else:
612 epsnew /= normv
613 F1=self.eval(x0 + epsnew * v)
614 return (F1-F0)/epsnew
615
616 ######################################
617 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):
618 """
619 solves a non-linear problem M{F(x)=0} for unknown M{x} using the stopping criterion:
620
621 M{norm(F(x) <= atol + rtol * norm(F(x0)}
622
623 where M{x0} is the initial guess.
624
625 @param defect: object defining the the function M{F}, C{defect.norm} defines the M{norm} used in the stopping criterion.
626 @type defect: L{Defect}
627 @param x: initial guess for the solution, C{x} is altered.
628 @type x: any object type allowing basic operations such as L{numarray.NumArray}, L{Data}
629 @param iter_max: maximum number of iteration steps
630 @type iter_max: positive C{int}
631 @param sub_iter_max:
632 @type sub_iter_max:
633 @param atol: absolute tolerance for the solution
634 @type atol: positive C{float}
635 @param rtol: relative tolerance for the solution
636 @type rtol: positive C{float}
637 @param gamma: tolerance safety factor for inner iteration
638 @type gamma: positive C{float}, less than 1
639 @param sub_tol_max: upper bound for inner tolerance.
640 @type sub_tol_max: positive C{float}, less than 1
641 @return: an approximation of the solution with the desired accuracy
642 @rtype: same type as the initial guess.
643 """
644 lmaxit=iter_max
645 if atol<0: raise ValueError,"atol needs to be non-negative."
646 if rtol<0: raise ValueError,"rtol needs to be non-negative."
647 if rtol+atol<=0: raise ValueError,"rtol or atol needs to be non-negative."
648 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
649 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
650
651 F=defect(x)
652 fnrm=defect.norm(F)
653 stop_tol=atol + rtol*fnrm
654 sub_tol=sub_tol_max
655 if verbose: print "NewtonGMRES: initial residual = %e."%fnrm
656 if verbose: print " tolerance = %e."%sub_tol
657 iter=1
658 #
659 # main iteration loop
660 #
661 while not fnrm<=stop_tol:
662 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
663 #
664 # adjust sub_tol_
665 #
666 if iter > 1:
667 rat=fnrm/fnrmo
668 sub_tol_old=sub_tol
669 sub_tol=gamma*rat**2
670 if gamma*sub_tol_old**2 > .1: sub_tol=max(sub_tol,gamma*sub_tol_old**2)
671 sub_tol=max(min(sub_tol,sub_tol_max), .5*stop_tol/fnrm)
672 #
673 # calculate newton increment xc
674 # if iter_max in __FDGMRES is reached MaxIterReached is thrown
675 # if iter_restart -1 is returned as sub_iter
676 # if atol is reached sub_iter returns the numer of steps performed to get there
677 #
678 #
679 if verbose: print " subiteration (GMRES) is called with relative tolerance %e."%sub_tol
680 try:
681 xc, sub_iter=__FDGMRES(F, defect, x, sub_tol*fnrm, iter_max=iter_max-iter, iter_restart=sub_iter_max)
682 except MaxIterReached:
683 raise MaxIterReached,"maximum number of %s steps reached."%iter_max
684 if sub_iter<0:
685 iter+=sub_iter_max
686 else:
687 iter+=sub_iter
688 # ====
689 x+=xc
690 F=defect(x)
691 iter+=1
692 fnrmo, fnrm=fnrm, defect.norm(F)
693 if verbose: print " step %s: residual %e."%(iter,fnrm)
694 if verbose: print "NewtonGMRES: completed after %s steps."%iter
695 return x
696
697 def __givapp(c,s,vin):
698 """
699 apply a sequence of Givens rotations (c,s) to the recuirsively to the vector vin
700 @warning: C{vin} is altered.
701 """
702 vrot=vin
703 if isinstance(c,float):
704 vrot=[c*vrot[0]-s*vrot[1],s*vrot[0]+c*vrot[1]]
705 else:
706 for i in range(len(c)):
707 w1=c[i]*vrot[i]-s[i]*vrot[i+1]
708 w2=s[i]*vrot[i]+c[i]*vrot[i+1]
709 vrot[i:i+2]=w1,w2
710 return vrot
711
712 def __FDGMRES(F0, defect, x0, atol, iter_max=100, iter_restart=20):
713 h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)
714 c=numarray.zeros(iter_restart,numarray.Float64)
715 s=numarray.zeros(iter_restart,numarray.Float64)
716 g=numarray.zeros(iter_restart,numarray.Float64)
717 v=[]
718
719 rho=defect.norm(F0)
720 if rho<=0.: return x0*0
721
722 v.append(-F0/rho)
723 g[0]=rho
724 iter=0
725 while rho > atol and iter<iter_restart-1:
726
727 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
728
729 p=defect.derivative(F0,x0,v[iter], v_is_normalised=True)
730 v.append(p)
731
732 v_norm1=defect.norm(v[iter+1])
733
734 # Modified Gram-Schmidt
735 for j in range(iter+1):
736 h[j,iter]=defect.bilinearform(v[j],v[iter+1])
737 v[iter+1]-=h[j,iter]*v[j]
738
739 h[iter+1,iter]=defect.norm(v[iter+1])
740 v_norm2=h[iter+1,iter]
741
742 # Reorthogonalize if needed
743 if v_norm1 + 0.001*v_norm2 == v_norm1: #Brown/Hindmarsh condition (default)
744 for j in range(iter+1):
745 hr=defect.bilinearform(v[j],v[iter+1])
746 h[j,iter]=h[j,iter]+hr
747 v[iter+1] -= hr*v[j]
748
749 v_norm2=defect.norm(v[iter+1])
750 h[iter+1,iter]=v_norm2
751 # watch out for happy breakdown
752 if not v_norm2 == 0:
753 v[iter+1]=v[iter+1]/h[iter+1,iter]
754
755 # Form and store the information for the new Givens rotation
756 if iter > 0 :
757 hhat=numarray.zeros(iter+1,numarray.Float64)
758 for i in range(iter+1) : hhat[i]=h[i,iter]
759 hhat=__givapp(c[0:iter],s[0:iter],hhat);
760 for i in range(iter+1) : h[i,iter]=hhat[i]
761
762 mu=math.sqrt(h[iter,iter]*h[iter,iter]+h[iter+1,iter]*h[iter+1,iter])
763
764 if mu!=0 :
765 c[iter]=h[iter,iter]/mu
766 s[iter]=-h[iter+1,iter]/mu
767 h[iter,iter]=c[iter]*h[iter,iter]-s[iter]*h[iter+1,iter]
768 h[iter+1,iter]=0.0
769 g[iter:iter+2]=__givapp(c[iter],s[iter],g[iter:iter+2])
770
771 # Update the residual norm
772 rho=abs(g[iter+1])
773 iter+=1
774
775 # At this point either iter > iter_max or rho < tol.
776 # It's time to compute x and leave.
777 if iter > 0 :
778 y=numarray.zeros(iter,numarray.Float64)
779 y[iter-1] = g[iter-1] / h[iter-1,iter-1]
780 if iter > 1 :
781 i=iter-2
782 while i>=0 :
783 y[i] = ( g[i] - numarray.dot(h[i,i+1:iter], y[i+1:iter])) / h[i,i]
784 i=i-1
785 xhat=v[iter-1]*y[iter-1]
786 for i in range(iter-1):
787 xhat += v[i]*y[i]
788 else :
789 xhat=v[0] * 0
790
791 if iter<iter_restart-1:
792 stopped=iter
793 else:
794 stopped=-1
795
796 return xhat,stopped
797
798 ##############################################
799 def GMRES(b, Aprod, bilinearform, atol=0, rtol=1.e-8, x=None, iter_max=100, iter_restart=20, verbose=False):
800 ################################################
801 m=iter_restart
802 iter=0
803 if x == None:
804 xc=b*0
805 else:
806 xc=x*1
807 if rtol>0:
808 r_dot_r = bilinearform(b, b)
809 if r_dot_r<0: raise NegativeNorm,"negative norm."
810 atol2=atol+rtol*math.sqrt(r_dot_r)
811 if verbose: print "GMRES: norm of right hand side = %e (absolute tolerance = %e)"%(math.sqrt(r_dot_r), atol2)
812 else:
813 atol2=atol
814 if verbose: print "GMRES: absolute tolerance = %e"%atol2
815
816 while True:
817 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached"%iter_max
818 xc,stopped=__GMRESm(b, Aprod, bilinearform, atol2, x=xc, iter_max=iter_max-iter, iter_restart=m, verbose=verbose)
819 iter+=iter_restart
820 if stopped: break
821 if verbose: print "GMRES: restart."
822 if verbose: print "GMRES: tolerance has reached."
823 return xc
824
825 def __GMRESm(b, Aprod, bilinearform, atol, x, iter_max=100, iter_restart=20, verbose=False):
826 iter=0
827
828 h=numarray.zeros((iter_restart+1,iter_restart),numarray.Float64)
829 c=numarray.zeros(iter_restart,numarray.Float64)
830 s=numarray.zeros(iter_restart,numarray.Float64)
831 g=numarray.zeros(iter_restart+1,numarray.Float64)
832 v=[]
833
834 r=b-Aprod(x)
835 r_dot_r = bilinearform(r, r)
836 if r_dot_r<0: raise NegativeNorm,"negative norm."
837 rho=math.sqrt(r_dot_r)
838
839 v.append(r/rho)
840 g[0]=rho
841
842 if verbose: print "GMRES: initial residual %e (absolute tolerance = %e)"%(rho,atol)
843 while not (rho<=atol or iter==iter_restart):
844
845 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
846
847 p=Aprod(v[iter])
848 v.append(p)
849
850 v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
851
852 # Modified Gram-Schmidt
853 for j in range(iter+1):
854 h[j,iter]=bilinearform(v[j],v[iter+1])
855 v[iter+1]-=h[j,iter]*v[j]
856
857 h[iter+1,iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
858 v_norm2=h[iter+1,iter]
859
860 # Reorthogonalize if needed
861 if v_norm1 + 0.001*v_norm2 == v_norm1: #Brown/Hindmarsh condition (default)
862 for j in range(iter+1):
863 hr=bilinearform(v[j],v[iter+1])
864 h[j,iter]=h[j,iter]+hr
865 v[iter+1] -= hr*v[j]
866
867 v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
868 h[iter+1,iter]=v_norm2
869
870 # watch out for happy breakdown
871 if not v_norm2 == 0:
872 v[iter+1]=v[iter+1]/h[iter+1,iter]
873
874 # Form and store the information for the new Givens rotation
875 if iter > 0: h[:iter+1,iter]=__givapp(c[:iter],s[:iter],h[:iter+1,iter])
876 mu=math.sqrt(h[iter,iter]*h[iter,iter]+h[iter+1,iter]*h[iter+1,iter])
877
878 if mu!=0 :
879 c[iter]=h[iter,iter]/mu
880 s[iter]=-h[iter+1,iter]/mu
881 h[iter,iter]=c[iter]*h[iter,iter]-s[iter]*h[iter+1,iter]
882 h[iter+1,iter]=0.0
883 g[iter:iter+2]=__givapp(c[iter],s[iter],g[iter:iter+2])
884 # Update the residual norm
885
886 rho=abs(g[iter+1])
887 if verbose: print "GMRES: iteration step %s: residual %e"%(iter,rho)
888 iter+=1
889
890 # At this point either iter > iter_max or rho < tol.
891 # It's time to compute x and leave.
892
893 if verbose: print "GMRES: iteration stopped after %s step."%iter
894 if iter > 0 :
895 y=numarray.zeros(iter,numarray.Float64)
896 y[iter-1] = g[iter-1] / h[iter-1,iter-1]
897 if iter > 1 :
898 i=iter-2
899 while i>=0 :
900 y[i] = ( g[i] - numarray.dot(h[i,i+1:iter], y[i+1:iter])) / h[i,i]
901 i=i-1
902 xhat=v[iter-1]*y[iter-1]
903 for i in range(iter-1):
904 xhat += v[i]*y[i]
905 else:
906 xhat=v[0] * 0
907
908 x += xhat
909 if iter<iter_restart-1:
910 stopped=True
911 else:
912 stopped=False
913
914 return x,stopped
915
916 #################################################
917 def MINRES(b, Aprod, Msolve, bilinearform, atol=0, rtol=1.e-8, x=None, iter_max=100):
918 #################################################
919 #
920 # minres solves the system of linear equations Ax = b
921 # where A is a symmetric matrix (possibly indefinite or singular)
922 # and b is a given vector.
923 #
924 # "A" may be a dense or sparse matrix (preferably sparse!)
925 # or the name of a function such that
926 # y = A(x)
927 # returns the product y = Ax for any given vector x.
928 #
929 # "M" defines a positive-definite preconditioner M = C C'.
930 # "M" may be a dense or sparse matrix (preferably sparse!)
931 # or the name of a function such that
932 # solves the system My = x for any given vector x.
933 #
934 #
935
936 #------------------------------------------------------------------
937 # Set up y and v for the first Lanczos vector v1.
938 # y = beta1 P' v1, where P = C**(-1).
939 # v is really P' v1.
940 #------------------------------------------------------------------
941 if x==None:
942 x=0*b
943 else:
944 b += (-1)*Aprod(x)
945
946 r1 = b
947 y = Msolve(b)
948 beta1 = bilinearform(y,b)
949
950 if beta1< 0: raise NegativeNorm,"negative norm."
951
952 # If b = 0 exactly, stop with x = 0.
953 if beta1==0: return x*0.
954
955 if beta1> 0:
956 beta1 = math.sqrt(beta1)
957
958 #------------------------------------------------------------------
959 # Initialize quantities.
960 # ------------------------------------------------------------------
961 iter = 0
962 Anorm = 0
963 ynorm = 0
964 oldb = 0
965 beta = beta1
966 dbar = 0
967 epsln = 0
968 phibar = beta1
969 rhs1 = beta1
970 rhs2 = 0
971 rnorm = phibar
972 tnorm2 = 0
973 ynorm2 = 0
974 cs = -1
975 sn = 0
976 w = b*0.
977 w2 = b*0.
978 r2 = r1
979 eps = 0.0001
980
981 #---------------------------------------------------------------------
982 # Main iteration loop.
983 # --------------------------------------------------------------------
984 while not rnorm<=atol+rtol*Anorm*ynorm: # checks ||r|| < (||A|| ||x||) * TOL
985
986 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
987 iter = iter + 1
988
989 #-----------------------------------------------------------------
990 # Obtain quantities for the next Lanczos vector vk+1, k = 1, 2,...
991 # The general iteration is similar to the case k = 1 with v0 = 0:
992 #
993 # p1 = Operator * v1 - beta1 * v0,
994 # alpha1 = v1'p1,
995 # q2 = p2 - alpha1 * v1,
996 # beta2^2 = q2'q2,
997 # v2 = (1/beta2) q2.
998 #
999 # Again, y = betak P vk, where P = C**(-1).
1000 #-----------------------------------------------------------------
1001 s = 1/beta # Normalize previous vector (in y).
1002 v = s*y # v = vk if P = I
1003
1004 y = Aprod(v)
1005
1006 if iter >= 2:
1007 y = y - (beta/oldb)*r1
1008
1009 alfa = bilinearform(v,y) # alphak
1010 y += (- alfa/beta)*r2
1011 r1 = r2
1012 r2 = y
1013 y = Msolve(r2)
1014 oldb = beta # oldb = betak
1015 beta = bilinearform(y,r2) # beta = betak+1^2
1016 if beta < 0: raise NegativeNorm,"negative norm."
1017
1018 beta = math.sqrt( beta )
1019 tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta
1020
1021 if iter==1: # Initialize a few things.
1022 gmax = abs( alfa ) # alpha1
1023 gmin = gmax # alpha1
1024
1025 # Apply previous rotation Qk-1 to get
1026 # [deltak epslnk+1] = [cs sn][dbark 0 ]
1027 # [gbar k dbar k+1] [sn -cs][alfak betak+1].
1028
1029 oldeps = epsln
1030 delta = cs * dbar + sn * alfa # delta1 = 0 deltak
1031 gbar = sn * dbar - cs * alfa # gbar 1 = alfa1 gbar k
1032 epsln = sn * beta # epsln2 = 0 epslnk+1
1033 dbar = - cs * beta # dbar 2 = beta2 dbar k+1
1034
1035 # Compute the next plane rotation Qk
1036
1037 gamma = math.sqrt(gbar*gbar+beta*beta) # gammak
1038 gamma = max(gamma,eps)
1039 cs = gbar / gamma # ck
1040 sn = beta / gamma # sk
1041 phi = cs * phibar # phik
1042 phibar = sn * phibar # phibark+1
1043
1044 # Update x.
1045
1046 denom = 1/gamma
1047 w1 = w2
1048 w2 = w
1049 w = (v - oldeps*w1 - delta*w2) * denom
1050 x += phi*w
1051
1052 # Go round again.
1053
1054 gmax = max(gmax,gamma)
1055 gmin = min(gmin,gamma)
1056 z = rhs1 / gamma
1057 ynorm2 = z*z + ynorm2
1058 rhs1 = rhs2 - delta*z
1059 rhs2 = - epsln*z
1060
1061 # Estimate various norms and test for convergence.
1062
1063 Anorm = math.sqrt( tnorm2 )
1064 ynorm = math.sqrt( ynorm2 )
1065
1066 rnorm = phibar
1067
1068 return x
1069
1070 def TFQMR(b, Aprod, Msolve, bilinearform, atol=0, rtol=1.e-8, x=None, iter_max=100):
1071
1072 # TFQMR solver for linear systems
1073 #
1074 #
1075 # initialization
1076 #
1077 errtol = math.sqrt(bilinearform(b,b))
1078 norm_b=errtol
1079 kmax = iter_max
1080 error = []
1081
1082 if math.sqrt(bilinearform(x,x)) != 0.0:
1083 r = b - Aprod(x)
1084 else:
1085 r = b
1086
1087 r=Msolve(r)
1088
1089 u1=0
1090 u2=0
1091 y1=0
1092 y2=0
1093
1094 w = r
1095 y1 = r
1096 iter = 0
1097 d = 0
1098
1099 v = Msolve(Aprod(y1))
1100 u1 = v
1101
1102 theta = 0.0;
1103 eta = 0.0;
1104 tau = math.sqrt(bilinearform(r,r))
1105 error = [ error, tau ]
1106 rho = tau * tau
1107 #
1108 # TFQMR iteration
1109 #
1110 # while ( iter < kmax-1 ):
1111
1112 while not tau<=atol+rtol*norm_b:
1113 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
1114
1115 sigma = bilinearform(r,v)
1116
1117 if ( sigma == 0.0 ):
1118 raise 'TFQMR breakdown, sigma=0'
1119
1120
1121 alpha = rho / sigma
1122
1123 for j in range(2):
1124 #
1125 # Compute y2 and u2 only if you have to
1126 #
1127 if ( j == 1 ):
1128 y2 = y1 - alpha * v
1129 u2 = Msolve(Aprod(y2))
1130
1131 m = 2 * (iter+1) - 2 + (j+1)
1132 if j==0:
1133 w = w - alpha * u1
1134 d = y1 + ( theta * theta * eta / alpha ) * d
1135 if j==1:
1136 w = w - alpha * u2
1137 d = y2 + ( theta * theta * eta / alpha ) * d
1138
1139 theta = math.sqrt(bilinearform(w,w))/ tau
1140 c = 1.0 / math.sqrt ( 1.0 + theta * theta )
1141 tau = tau * theta * c
1142 eta = c * c * alpha
1143 x = x + eta * d
1144 #
1145 # Try to terminate the iteration at each pass through the loop
1146 #
1147 # if ( tau * math.sqrt ( m + 1 ) <= errtol ):
1148 # error = [ error, tau ]
1149 # total_iters = iter
1150 # break
1151
1152
1153 if ( rho == 0.0 ):
1154 raise 'TFQMR breakdown, rho=0'
1155
1156
1157 rhon = bilinearform(r,w)
1158 beta = rhon / rho;
1159 rho = rhon;
1160 y1 = w + beta * y2;
1161 u1 = Msolve(Aprod(y1))
1162 v = u1 + beta * ( u2 + beta * v )
1163 error = [ error, tau ]
1164 total_iters = iter
1165
1166 iter = iter + 1
1167
1168 return x
1169
1170
1171 #############################################
1172
1173 class ArithmeticTuple(object):
1174 """
1175 tuple supporting inplace update x+=y and scaling x=a*y where x,y is an ArithmeticTuple and a is a float.
1176
1177 example of usage:
1178
1179 from esys.escript import Data
1180 from numarray import array
1181 a=Data(...)
1182 b=array([1.,4.])
1183 x=ArithmeticTuple(a,b)
1184 y=5.*x
1185
1186 """
1187 def __init__(self,*args):
1188 """
1189 initialize object with elements args.
1190
1191 @param args: tuple of object that support implace add (x+=y) and scaling (x=a*y)
1192 """
1193 self.__items=list(args)
1194
1195 def __len__(self):
1196 """
1197 number of items
1198
1199 @return: number of items
1200 @rtype: C{int}
1201 """
1202 return len(self.__items)
1203
1204 def __getitem__(self,index):
1205 """
1206 get an item
1207
1208 @param index: item to be returned
1209 @type index: C{int}
1210 @return: item with index C{index}
1211 """
1212 return self.__items.__getitem__(index)
1213
1214 def __mul__(self,other):
1215 """
1216 scaling from the right
1217
1218 @param other: scaling factor
1219 @type other: C{float}
1220 @return: itemwise self*other
1221 @rtype: L{ArithmeticTuple}
1222 """
1223 out=[]
1224 try:
1225 l=len(other)
1226 if l!=len(self):
1227 raise ValueError,"length of of arguments don't match."
1228 for i in range(l): out.append(self[i]*other[i])
1229 except TypeError:
1230 for i in range(len(self)): out.append(self[i]*other)
1231 return ArithmeticTuple(*tuple(out))
1232
1233 def __rmul__(self,other):
1234 """
1235 scaling from the left
1236
1237 @param other: scaling factor
1238 @type other: C{float}
1239 @return: itemwise other*self
1240 @rtype: L{ArithmeticTuple}
1241 """
1242 out=[]
1243 try:
1244 l=len(other)
1245 if l!=len(self):
1246 raise ValueError,"length of of arguments don't match."
1247 for i in range(l): out.append(other[i]*self[i])
1248 except TypeError:
1249 for i in range(len(self)): out.append(other*self[i])
1250 return ArithmeticTuple(*tuple(out))
1251
1252 def __div__(self,other):
1253 """
1254 dividing from the right
1255
1256 @param other: scaling factor
1257 @type other: C{float}
1258 @return: itemwise self/other
1259 @rtype: L{ArithmeticTuple}
1260 """
1261 return self*(1/other)
1262
1263 def __rdiv__(self,other):
1264 """
1265 dividing from the left
1266
1267 @param other: scaling factor
1268 @type other: C{float}
1269 @return: itemwise other/self
1270 @rtype: L{ArithmeticTuple}
1271 """
1272 out=[]
1273 try:
1274 l=len(other)
1275 if l!=len(self):
1276 raise ValueError,"length of of arguments don't match."
1277 for i in range(l): out.append(other[i]/self[i])
1278 except TypeError:
1279 for i in range(len(self)): out.append(other/self[i])
1280 return ArithmeticTuple(*tuple(out))
1281
1282 def __iadd__(self,other):
1283 """
1284 in-place add of other to self
1285
1286 @param other: increment
1287 @type other: C{ArithmeticTuple}
1288 """
1289 if len(self) != len(other):
1290 raise ValueError,"tuple length must match."
1291 for i in range(len(self)):
1292 self.__items[i]+=other[i]
1293 return self
1294
1295 def __add__(self,other):
1296 """
1297 add other to self
1298
1299 @param other: increment
1300 @type other: C{ArithmeticTuple}
1301 """
1302 out=[]
1303 try:
1304 l=len(other)
1305 if l!=len(self):
1306 raise ValueError,"length of of arguments don't match."
1307 for i in range(l): out.append(self[i]+other[i])
1308 except TypeError:
1309 for i in range(len(self)): out.append(self[i]+other)
1310 return ArithmeticTuple(*tuple(out))
1311
1312 def __sub__(self,other):
1313 """
1314 subtract other from self
1315
1316 @param other: increment
1317 @type other: C{ArithmeticTuple}
1318 """
1319 out=[]
1320 try:
1321 l=len(other)
1322 if l!=len(self):
1323 raise ValueError,"length of of arguments don't match."
1324 for i in range(l): out.append(self[i]-other[i])
1325 except TypeError:
1326 for i in range(len(self)): out.append(self[i]-other)
1327 return ArithmeticTuple(*tuple(out))
1328
1329 def __isub__(self,other):
1330 """
1331 subtract other from self
1332
1333 @param other: increment
1334 @type other: C{ArithmeticTuple}
1335 """
1336 if len(self) != len(other):
1337 raise ValueError,"tuple length must match."
1338 for i in range(len(self)):
1339 self.__items[i]-=other[i]
1340 return self
1341
1342 def __neg__(self):
1343 """
1344 negate
1345
1346 """
1347 out=[]
1348 for i in range(len(self)):
1349 out.append(-self[i])
1350 return ArithmeticTuple(*tuple(out))
1351
1352
1353 class HomogeneousSaddlePointProblem(object):
1354 """
1355 This provides a framwork for solving linear homogeneous saddle point problem of the form
1356
1357 Av+B^*p=f
1358 Bv =0
1359
1360 for the unknowns v and p and given operators A and B and given right hand side f.
1361 B^* is the adjoint operator of B.
1362 """
1363 def __init__(self,**kwargs):
1364 self.setTolerance()
1365 self.setAbsoluteTolerance()
1366 self.setSubToleranceReductionFactor()
1367
1368 #=============================================================
1369 def initialize(self):
1370 """
1371 initialize the problem (overwrite)
1372 """
1373 pass
1374
1375 def B(self,v):
1376 """
1377 returns Bv (overwrite)
1378 @rtype: equal to the type of p
1379
1380 @note: boundary conditions on p should be zero!
1381 """
1382 raise NotImplementedError,"no operator B implemented."
1383
1384 def inner_pBv(self,p,Bv):
1385 """
1386 returns inner product of element p and Bv (overwrite)
1387
1388 @type p: equal to the type of p
1389 @type Bv: equal to the type of result of operator B
1390 @rtype: C{float}
1391
1392 @rtype: equal to the type of p
1393 """
1394 raise NotImplementedError,"no inner product for p implemented."
1395
1396 def inner_p(self,p0,p1):
1397 """
1398 returns inner product of element p0 and p1 (overwrite)
1399
1400 @type p0: equal to the type of p
1401 @type p1: equal to the type of p
1402 @rtype: equal to the type of p
1403 """
1404 raise NotImplementedError,"no inner product for p implemented."
1405
1406 def inner_v(self,v0,v1):
1407 """
1408 returns inner product of two element v0 and v1 (overwrite)
1409
1410 @type v0: equal to the type of v
1411 @type v1: equal to the type of v
1412 @rtype: C{float}
1413
1414 @rtype: equal to the type of v
1415 """
1416 raise NotImplementedError,"no inner product for v implemented."
1417 pass
1418
1419 def solve_A(self,u,p):
1420 """
1421 solves Av=f-Au-B^*p with accuracy self.getSubProblemTolerance() (overwrite)
1422
1423 @rtype: equal to the type of v
1424 @note: boundary conditions on v should be zero!
1425 """
1426 raise NotImplementedError,"no operator A implemented."
1427
1428 def solve_prec(self,p):
1429 """
1430 provides a preconditioner for BA^{-1}B^* with accuracy self.getSubProblemTolerance() (overwrite)
1431
1432 @rtype: equal to the type of p
1433 @note: boundary conditions on p should be zero!
1434 """
1435 raise NotImplementedError,"no preconditioner for Schur complement implemented."
1436 #=============================================================
1437 def __Aprod_PCG(self,p):
1438 # return (v,Bv) with v=A^-1B*p
1439 #solve Av =B^*p as Av =f-Az-B^*(-p)
1440 v=self.solve_A(self.__z,-p)
1441 return ArithmeticTuple(v, self.B(v))
1442
1443 def __inner_PCG(self,p,r):
1444 return self.inner_pBv(p,r[1])
1445
1446 def __Msolve_PCG(self,r):
1447 return self.solve_prec(r[1])
1448 #=============================================================
1449 def __Aprod_GMRES(self,x):
1450 # return w,q from v, p
1451 # solve Aw =Av+B^*p as Aw =f-A(z-v)-B^*(-p)
1452 # and Sq=B(v-w)
1453 v=x[0]
1454 p=x[1]
1455 w=self.solve_A(self.__z-v,-p)
1456 Bw=self.B(w-v)
1457 q=self.solve_prec(Bw)
1458 return ArithmeticTuple(w,q)
1459
1460 def __inner_GMRES(self,a1,a2):
1461 return self.inner_p(a1[1],a2[1])+self.inner_v(a1[0],a2[0])
1462
1463 #=============================================================
1464 def norm(self,v,p):
1465 f=self.inner_p(p,p)+self.inner_v(v,v)
1466 if f<0:
1467 raise ValueError,"negative norm."
1468 return math.sqrt(f)
1469
1470 def solve(self,v,p,max_iter=20, verbose=False, show_details=False, useUzawa=True, iter_restart=20,max_correction_steps=3):
1471 """
1472 solves the saddle point problem using initial guesses v and p.
1473
1474 @param v: initial guess for velocity
1475 @param p: initial guess for pressure
1476 @type v: L{Data}
1477 @type p: L{Data}
1478 @param useUzawa: indicate the usage of the Uzawa scheme. Otherwise the problem is solved in coupled form.
1479 @param max_iter: maximum number of iteration steps per correction attempt.
1480 @param verbose: show information on the progress of the saddlepoint problem solver.
1481 @param show_details: show details of the sub problems solves
1482 @param iter_restart: restart the iteration after C{iter_restart} steps (only used if useUzaw=False)
1483 @param max_correction_steps: maximum number of iteration steps in the attempt get |Bv| to zero.
1484 @return: new approximations for velocity and pressure
1485 @type useUzawa: C{bool}
1486 @type max_iter: C{int}
1487 @type verbose: C{bool}
1488 @type show_details: C{bool}
1489 @type iter_restart: C{int}
1490 @type max_correction_steps: C{int}
1491 @rtype: C{tuple} of L{Data} objects
1492 """
1493 self.verbose=verbose
1494 self.show_details=show_details and self.verbose
1495 #=================================================================================
1496 # Az=f is solved as A(z-v)=f-Av-B^*0 (typically with z-v=0 on boundary)
1497 self.setSubProblemTolerance(self.getTolerance(), apply_reduction=True)
1498 self.__z=v+self.solve_A(v,p*0)
1499 Bz=self.B(self.__z)
1500 p0=self.solve_prec(Bz)
1501 # tolerances:
1502 rtol=self.getTolerance()
1503 atol=self.getAbsoluteTolerance()
1504 f0=self.norm(self.__z,p0)
1505 if not f0>0:
1506 return self.__z, p*0
1507 ATOL=rtol*f0+atol
1508 if not ATOL>0:
1509 raise ValueError,"overall absolute tolerance needs to be positive."
1510 if self.verbose: print "saddle point solver: initial residual %e, tolerance = %e."%(f0,ATOL)
1511 # initialization
1512 self.iter=0
1513 correction_step=0
1514 converged=False
1515 # initial guess:
1516 q=p*1
1517 if not useUzawa:
1518 u=v
1519 while not converged :
1520 if useUzawa:
1521 self.setSubProblemTolerance(ATOL/f0, apply_reduction=True)
1522 # assume p is known: then v=z-A^-1B^*p
1523 #
1524 # which leads to BA^-1B^*p = Bz
1525 #
1526 # note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bu
1527 #
1528 # note that q=BA^-1B^*p = B (A^-1B^*p)=Bw with Aw=B^*p = f - Az - B^*(-p)
1529 #
1530 if self.verbose: print "enter PCG method (iter_max=%s, atol=%e, subtolerance=%e)"%(max_iter,ATOL, self.getSubProblemTolerance())
1531 q,r=PCG(ArithmeticTuple(self.__z,Bz),self.__Aprod_PCG,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL, rtol=0.,iter_max=max_iter, x=q, verbose=self.verbose)
1532 u=r[0]
1533 Bu=r[1]
1534 else:
1535 self.setSubProblemTolerance(ATOL/f0, apply_reduction=True)
1536 #
1537 # with v=dv+z
1538 #
1539 # A dv + B p = A (v-z) + Bp = 0
1540 # B dv = - Bz
1541 #
1542 # apply the preconditioner [[A^{-1} 0][(S^{-1} B A^{-1}) -S^{-1}]] to the problem
1543 # the right hand side is then [0, S^{-1} Bz ]
1544 #
1545 #
1546 if self.verbose: print "enter GMRES (iter_max=%s, atol=%e, subtolerance=%e)"%(max_iter,ATOL,self.getSubProblemTolerance())
1547 x=GMRES(ArithmeticTuple(0*self.__z,p0),self.__Aprod_GMRES,self.__inner_GMRES,atol=ATOL, rtol=0.,iter_max=max_iter, x=ArithmeticTuple(u-self.__z,q),iter_restart=iter_restart, verbose=self.verbose)
1548 self.setSubProblemTolerance(ATOL/f0, apply_reduction=False)
1549 u=self.__z+x[0]
1550 q=x[1]
1551 Bu=self.B(u)
1552 # now we check if |Bu| ~ 0 or more precise |Bu| <= rtol * |v|
1553 nu=self.inner_v(u,u)
1554 nz=self.inner_v(self.__z,self.__z)
1555 p2=self.solve_prec(Bu)
1556 nBu=self.inner_p(p2,p2)
1557 if not nu>=0 and not nBu>=0: raise NegativeNorm,"negative norm."
1558 nu=math.sqrt(nu)
1559 nBu=math.sqrt(nBu)
1560 if self.verbose: print "saddle point solver: norm v= %e (Bv = %e)"%(nu,nBu)
1561 QTOL=atol+nu*rtol
1562 if nBu <= QTOL:
1563 converged=True
1564 else:
1565 ATOL=QTOL/nBu*ATOL*0.3
1566 if self.verbose: print "correction step %s: tolerance reduced to %e."%(correction_step,ATOL)
1567 converged=False
1568 correction_step+=1
1569 if correction_step>max_correction_steps:
1570 raise CorrectionFailed,"Given up after %d correction steps."%correction_step
1571 if self.verbose: print "saddle point solver: tolerance reached."
1572 return u,q
1573
1574 #==========================================================================================================================
1575 def setTolerance(self,tolerance=1.e-4):
1576 """
1577 sets the relative tolerance for (v,p)
1578
1579 @param tolerance: tolerance to be used
1580 @type tolerance: non-negative C{float}
1581 """
1582 if tolerance<0:
1583 raise ValueError,"tolerance must be positive."
1584 self.__rtol=tolerance
1585 def getTolerance(self):
1586 """
1587 returns the relative tolerance
1588
1589 @return: relative tolerance
1590 @rtype: C{float}
1591 """
1592 return self.__rtol
1593 def setAbsoluteTolerance(self,tolerance=0.):
1594 """
1595 sets the absolute tolerance
1596
1597 @param tolerance: tolerance to be used
1598 @type tolerance: non-negative C{float}
1599 """
1600 if tolerance<0:
1601 raise ValueError,"tolerance must be non-negative."
1602 self.__atol=tolerance
1603 def getAbsoluteTolerance(self):
1604 """
1605 returns the absolute tolerance
1606
1607 @return: absolute tolerance
1608 @rtype: C{float}
1609 """
1610 return self.__atol
1611 def setSubToleranceReductionFactor(self,reduction=None):
1612 """
1613 sets the reduction factor for the tolerance used to solve the subproblem in each iteration step.
1614 if set to C{None} the square of the problem tolerance is used.
1615
1616 @param reduction: reduction factor.
1617 @type reduction: positive C{float} less or equal one or C{None}
1618 """
1619 if not reduction == None:
1620 if reduction<=0 or reduction>1:
1621 raise ValueError,"reduction factor (=%e) must be positive and less or equal one."%reduction
1622 self.__reduction=reduction
1623 def getSubToleranceReductionFactor(self):
1624 """
1625 returns the subproblem reduction factor
1626
1627 @return: subproblem reduction factor
1628 @rtype: C{float} or C{None}
1629 """
1630 return self.__reduction
1631
1632 def setSubProblemTolerance(self,rtol,apply_reduction=False):
1633 """
1634 sets the relative tolerance to solve the subproblem(s).
1635
1636 @param rtol: relative tolerence
1637 @type rtol: positive C{float}
1638 @param apply_reduction: if set the the reduction factor is applied.
1639 @type apply_reduction: C{bool}
1640 """
1641 if rtol<=0:
1642 raise ValueError,"tolerance must be positive."
1643 if apply_reduction:
1644 if self.__reduction == None:
1645 self.__sub_tol=max(rtol**2,util.EPSILON*10)
1646 else:
1647 self.__sub_tol=max(self.__reduction*rtol,util.EPSILON*10)
1648 else:
1649 self.__sub_tol=rtol
1650 def getSubProblemTolerance(self):
1651 """
1652 returns the subproblem reduction factor
1653
1654 @return: subproblem reduction factor
1655 @rtype: C{float}
1656 """
1657 return self.__sub_tol
1658
1659 def MaskFromBoundaryTag(domain,*tags):
1660 """
1661 create a mask on the Solution(domain) function space which one for samples
1662 that touch the boundary tagged by tags.
1663
1664 usage: m=MaskFromBoundaryTag(domain,"left", "right")
1665
1666 @param domain: a given domain
1667 @type domain: L{escript.Domain}
1668 @param tags: boundray tags
1669 @type tags: C{str}
1670 @return: a mask which marks samples that are touching the boundary tagged by any of the given tags.
1671 @rtype: L{escript.Data} of rank 0
1672 """
1673 pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1674 d=escript.Scalar(0.,escript.FunctionOnBoundary(domain))
1675 for t in tags: d.setTaggedValue(t,1.)
1676 pde.setValue(y=d)
1677 return util.whereNonZero(pde.getRightHandSide())
1678 #============================================================================================================================================
1679 class SaddlePointProblem(object):
1680 """
1681 This implements a solver for a saddlepoint problem
1682
1683 M{f(u,p)=0}
1684 M{g(u)=0}
1685
1686 for u and p. The problem is solved with an inexact Uszawa scheme for p:
1687
1688 M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})}
1689 M{Q_g (p^{k+1}-p^{k}) = g(u^{k+1})}
1690
1691 where Q_f is an approximation of the Jacobiean A_f of f with respect to u and Q_f is an approximation of
1692 A_g A_f^{-1} A_g with A_g is the jacobiean of g with respect to p. As a the construction of a 'proper'
1693 Q_g can be difficult, non-linear conjugate gradient method is applied to solve for p, so Q_g plays
1694 in fact the role of a preconditioner.
1695 """
1696 def __init__(self,verbose=False,*args):
1697 """
1698 initializes the problem
1699
1700 @param verbose: switches on the printing out some information
1701 @type verbose: C{bool}
1702 @note: this method may be overwritten by a particular saddle point problem
1703 """
1704 print "SaddlePointProblem should not be used anymore!"
1705 if not isinstance(verbose,bool):
1706 raise TypeError("verbose needs to be of type bool.")
1707 self.__verbose=verbose
1708 self.relaxation=1.
1709 DeprecationWarning("SaddlePointProblem should not be used anymore.")
1710
1711 def trace(self,text):
1712 """
1713 prints text if verbose has been set
1714
1715 @param text: a text message
1716 @type text: C{str}
1717 """
1718 if self.__verbose: print "%s: %s"%(str(self),text)
1719
1720 def solve_f(self,u,p,tol=1.e-8):
1721 """
1722 solves
1723
1724 A_f du = f(u,p)
1725
1726 with tolerance C{tol} and return du. A_f is Jacobiean of f with respect to u.
1727
1728 @param u: current approximation of u
1729 @type u: L{escript.Data}
1730 @param p: current approximation of p
1731 @type p: L{escript.Data}
1732 @param tol: tolerance expected for du
1733 @type tol: C{float}
1734 @return: increment du
1735 @rtype: L{escript.Data}
1736 @note: this method has to be overwritten by a particular saddle point problem
1737 """
1738 pass
1739
1740 def solve_g(self,u,tol=1.e-8):
1741 """
1742 solves
1743
1744 Q_g dp = g(u)
1745
1746 with Q_g is a preconditioner for A_g A_f^{-1} A_g with A_g is the jacobiean of g with respect to p.
1747
1748 @param u: current approximation of u
1749 @type u: L{escript.Data}
1750 @param tol: tolerance expected for dp
1751 @type tol: C{float}
1752 @return: increment dp
1753 @rtype: L{escript.Data}
1754 @note: this method has to be overwritten by a particular saddle point problem
1755 """
1756 pass
1757
1758 def inner(self,p0,p1):
1759 """
1760 inner product of p0 and p1 approximating p. Typically this returns integrate(p0*p1)
1761 @return: inner product of p0 and p1
1762 @rtype: C{float}
1763 """
1764 pass
1765
1766 subiter_max=3
1767 def solve(self,u0,p0,tolerance=1.e-6,tolerance_u=None,iter_max=100,accepted_reduction=0.995,relaxation=None):
1768 """
1769 runs the solver
1770
1771 @param u0: initial guess for C{u}
1772 @type u0: L{esys.escript.Data}
1773 @param p0: initial guess for C{p}
1774 @type p0: L{esys.escript.Data}
1775 @param tolerance: tolerance for relative error in C{u} and C{p}
1776 @type tolerance: positive C{float}
1777 @param tolerance_u: tolerance for relative error in C{u} if different from C{tolerance}
1778 @type tolerance_u: positive C{float}
1779 @param iter_max: maximum number of iteration steps.
1780 @type iter_max: C{int}
1781 @param accepted_reduction: if the norm g cannot be reduced by C{accepted_reduction} backtracking to adapt the
1782 relaxation factor. If C{accepted_reduction=None} no backtracking is used.
1783 @type accepted_reduction: positive C{float} or C{None}
1784 @param relaxation: initial relaxation factor. If C{relaxation==None}, the last relaxation factor is used.
1785 @type relaxation: C{float} or C{None}
1786 """
1787 tol=1.e-2
1788 if tolerance_u==None: tolerance_u=tolerance
1789 if not relaxation==None: self.relaxation=relaxation
1790 if accepted_reduction ==None:
1791 angle_limit=0.
1792 elif accepted_reduction>=1.:
1793 angle_limit=0.
1794 else:
1795 angle_limit=util.sqrt(1-accepted_reduction**2)
1796 self.iter=0
1797 u=u0
1798 p=p0
1799 #
1800 # initialize things:
1801 #
1802 converged=False
1803 #
1804 # start loop:
1805 #
1806 # initial search direction is g
1807 #
1808 while not converged :
1809 if self.iter>iter_max:
1810 raise ArithmeticError("no convergence after %s steps."%self.iter)
1811 f_new=self.solve_f(u,p,tol)
1812 norm_f_new = util.Lsup(f_new)
1813 u_new=u-f_new
1814 g_new=self.solve_g(u_new,tol)
1815 self.iter+=1
1816 norm_g_new = util.sqrt(self.inner(g_new,g_new))
1817 if norm_f_new==0. and norm_g_new==0.: return u, p
1818 if self.iter>1 and not accepted_reduction==None:
1819 #
1820 # did we manage to reduce the norm of G? I
1821 # if not we start a backtracking procedure
1822 #
1823 # print "new/old norm = ",norm_g_new, norm_g, norm_g_new/norm_g
1824 if norm_g_new > accepted_reduction * norm_g:
1825 sub_iter=0
1826 s=self.relaxation
1827 d=g
1828 g_last=g
1829 self.trace(" start substepping: f = %s, g = %s, relaxation = %s."%(norm_f_new, norm_g_new, s))
1830 while sub_iter < self.subiter_max and norm_g_new > accepted_reduction * norm_g:
1831 dg= g_new-g_last
1832 norm_dg=abs(util.sqrt(self.inner(dg,dg))/self.relaxation)
1833 rad=self.inner(g_new,dg)/self.relaxation
1834 # print " ",sub_iter,": rad, norm_dg:",abs(rad), norm_dg*norm_g_new * angle_limit
1835 # print " ",sub_iter,": rad, norm_dg:",rad, norm_dg, norm_g_new, norm_g
1836 if abs(rad) < norm_dg*norm_g_new * angle_limit:
1837 if sub_iter>0: self.trace(" no further improvements expected from backtracking.")
1838 break
1839 r=self.relaxation
1840 self.relaxation= - rad/norm_dg**2
1841 s+=self.relaxation
1842 #####
1843 # a=g_new+self.relaxation*dg/r
1844 # print "predicted new norm = ",util.sqrt(self.inner(a,a)),util.sqrt(self.inner(g_new,g_new)), self.relaxation
1845 #####
1846 g_last=g_new
1847 p+=self.relaxation*d
1848 f_new=self.solve_f(u,p,tol)
1849 u_new=u-f_new
1850 g_new=self.solve_g(u_new,tol)
1851 self.iter+=1
1852 norm_f_new = util.Lsup(f_new)
1853 norm_g_new = util.sqrt(self.inner(g_new,g_new))
1854 # print " ",sub_iter," new g norm",norm_g_new
1855 self.trace(" %s th sub-step: f = %s, g = %s, relaxation = %s."%(sub_iter, norm_f_new, norm_g_new, s))
1856 #
1857 # can we expect reduction of g?
1858 #
1859 # u_last=u_new
1860 sub_iter+=1
1861 self.relaxation=s
1862 #
1863 # check for convergence:
1864 #
1865 norm_u_new = util.Lsup(u_new)
1866 p_new=p+self.relaxation*g_new
1867 norm_p_new = util.sqrt(self.inner(p_new,p_new))
1868 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))
1869
1870 if self.iter>1:
1871 dg2=g_new-g
1872 df2=f_new-f
1873 norm_dg2=util.sqrt(self.inner(dg2,dg2))
1874 norm_df2=util.Lsup(df2)
1875 # print norm_g_new, norm_g, norm_dg, norm_p, tolerance
1876 tol_eq_g=tolerance*norm_dg2/(norm_g*abs(self.relaxation))*norm_p_new
1877 tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new
1878 if norm_g_new <= tol_eq_g and norm_f_new <= tol_eq_f:
1879 converged=True
1880 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
1881 self.trace("convergence after %s steps."%self.iter)
1882 return u,p

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26