/[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 1896 - (show annotations)
Sun Oct 19 23:16:21 2008 UTC (10 years, 11 months ago) by gross
File MIME type: text/x-python
File size: 59627 byte(s)
small fixes.
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 IterationBreakDown(SolverSchemeException):
430 """
431 iteration scheme econouters an incurable breakdown.
432 """
433 pass
434 class NegativeNorm(SolverSchemeException):
435 """
436 a norm calculation returns a negative norm.
437 """
438 pass
439
440 class IterationHistory(object):
441 """
442 The IterationHistory class is used to define a stopping criterium. It keeps track of the
443 residual norms. The stoppingcriterium indicates termination if the residual norm has been reduced by
444 a given tolerance.
445 """
446 def __init__(self,tolerance=math.sqrt(util.EPSILON),verbose=False):
447 """
448 Initialization
449
450 @param tolerance: tolerance
451 @type tolerance: positive C{float}
452 @param verbose: switches on the printing out some information
453 @type verbose: C{bool}
454 """
455 if not tolerance>0.:
456 raise ValueError,"tolerance needs to be positive."
457 self.tolerance=tolerance
458 self.verbose=verbose
459 self.history=[]
460 def stoppingcriterium(self,norm_r,r,x):
461 """
462 returns True if the C{norm_r} is C{tolerance}*C{norm_r[0]} where C{norm_r[0]} is the residual norm at the first call.
463
464
465 @param norm_r: current residual norm
466 @type norm_r: non-negative C{float}
467 @param r: current residual (not used)
468 @param x: current solution approximation (not used)
469 @return: C{True} is the stopping criterium is fullfilled. Otherwise C{False} is returned.
470 @rtype: C{bool}
471
472 """
473 self.history.append(norm_r)
474 if self.verbose: print "iter: %s: inner(rhat,r) = %e"#(len(self.history)-1, self.history[-1])
475 return self.history[-1]<=self.tolerance * self.history[0]
476
477 def stoppingcriterium2(self,norm_r,norm_b,solver="GMRES",TOL=None):
478 """
479 returns True if the C{norm_r} is C{tolerance}*C{norm_b}
480
481
482 @param norm_r: current residual norm
483 @type norm_r: non-negative C{float}
484 @param norm_b: norm of right hand side
485 @type norm_b: non-negative C{float}
486 @return: C{True} is the stopping criterium is fullfilled. Otherwise C{False} is returned.
487 @rtype: C{bool}
488
489 """
490 if TOL==None:
491 TOL=self.tolerance
492 self.history.append(norm_r)
493 if self.verbose: print "iter: %s: norm(r) = %e"#(len(self.history)-1, self.history[-1])
494 return self.history[-1]<=TOL * norm_b
495
496 def PCG(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
497 """
498 Solver for
499
500 M{Ax=b}
501
502 with a symmetric and positive definite operator A (more details required!).
503 It uses the conjugate gradient method with preconditioner M providing an approximation of A.
504
505 The iteration is terminated if the C{stoppingcriterium} function return C{True}.
506
507 For details on the preconditioned conjugate gradient method see the book:
508
509 Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
510 T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
511 C. Romine, and H. van der Vorst.
512
513 @param b: the right hand side of the liner system. C{b} is altered.
514 @type b: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
515 @param Aprod: returns the value Ax
516 @type Aprod: function C{Aprod(x)} where C{x} is of the same object like argument C{x}. The returned object needs to be of the same type like argument C{b}.
517 @param Msolve: solves Mx=r
518 @type Msolve: function C{Msolve(r)} where C{r} is of the same type like argument C{b}. The returned object needs to be of the same
519 type like argument C{x}.
520 @param bilinearform: inner product C{<x,r>}
521 @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same type like argument C{x} and C{r} is . The returned value is a C{float}.
522 @param stoppingcriterium: function which returns True if a stopping criterium is meet. C{stoppingcriterium} has the arguments C{norm_r}, C{r} and C{x} giving the current norm of the residual (=C{sqrt(bilinearform(Msolve(r),r)}), the current residual and the current solution approximation. C{stoppingcriterium} is called in each iteration step.
523 @type stoppingcriterium: function that returns C{True} or C{False}
524 @param x: an initial guess for the solution. If no C{x} is given 0*b is used.
525 @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
526 @param iter_max: maximum number of iteration steps.
527 @type iter_max: C{int}
528 @return: the solution approximation and the corresponding residual
529 @rtype: C{tuple}
530 @warning: C{b} and C{x} are altered.
531 """
532 iter=0
533 if x==None:
534 x=0*b
535 else:
536 b += (-1)*Aprod(x)
537 r=b
538 rhat=Msolve(r)
539 d = rhat
540 rhat_dot_r = bilinearform(rhat, r)
541 if rhat_dot_r<0: raise NegativeNorm,"negative norm."
542
543 while not stoppingcriterium(math.sqrt(rhat_dot_r),r,x):
544 iter+=1
545 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
546
547 q=Aprod(d)
548 alpha = rhat_dot_r / bilinearform(d, q)
549 x += alpha * d
550 r += (-alpha) * q
551
552 rhat=Msolve(r)
553 rhat_dot_r_new = bilinearform(rhat, r)
554 beta = rhat_dot_r_new / rhat_dot_r
555 rhat+=beta * d
556 d=rhat
557
558 rhat_dot_r = rhat_dot_r_new
559 if rhat_dot_r<0: raise NegativeNorm,"negative norm."
560
561 return x,r
562
563 class Defect(object):
564 """
565 defines a non-linear defect F(x) of a variable x
566 """
567 def __init__(self):
568 """
569 initialize defect
570 """
571 self.setDerivativeIncrementLength()
572
573 def bilinearform(self, x0, x1):
574 """
575 returns the inner product of x0 and x1
576 @param x0: a value for x
577 @param x1: a value for x
578 @return: the inner product of x0 and x1
579 @rtype: C{float}
580 """
581 return 0
582
583 def norm(self,x):
584 """
585 the norm of argument C{x}
586
587 @param x: a value for x
588 @return: norm of argument x
589 @rtype: C{float}
590 @note: by default C{sqrt(self.bilinearform(x,x)} is retrurned.
591 """
592 s=self.bilinearform(x,x)
593 if s<0: raise NegativeNorm,"negative norm."
594 return math.sqrt(s)
595
596
597 def eval(self,x):
598 """
599 returns the value F of a given x
600
601 @param x: value for which the defect C{F} is evalulated.
602 @return: value of the defect at C{x}
603 """
604 return 0
605
606 def __call__(self,x):
607 return self.eval(x)
608
609 def setDerivativeIncrementLength(self,inc=math.sqrt(util.EPSILON)):
610 """
611 sets the relative length of the increment used to approximate the derivative of the defect
612 the increment is inc*norm(x)/norm(v)*v in the direction of v with x as a staring point.
613
614 @param inc: relative increment length
615 @type inc: positive C{float}
616 """
617 if inc<=0: raise ValueError,"positive increment required."
618 self.__inc=inc
619
620 def getDerivativeIncrementLength(self):
621 """
622 returns the relative increment length used to approximate the derivative of the defect
623 @return: value of the defect at C{x}
624 @rtype: positive C{float}
625 """
626 return self.__inc
627
628 def derivative(self, F0, x0, v, v_is_normalised=True):
629 """
630 returns the directional derivative at x0 in the direction of v
631
632 @param F0: value of this defect at x0
633 @param x0: value at which derivative is calculated.
634 @param v: direction
635 @param v_is_normalised: is true to indicate that C{v} is nomalized (self.norm(v)=0)
636 @return: derivative of this defect at x0 in the direction of C{v}
637 @note: by default numerical evaluation (self.eval(x0+eps*v)-F0)/eps is used but this method
638 maybe oepsnew verwritten to use exact evalution.
639 """
640 normx=self.norm(x0)
641 if normx>0:
642 epsnew = self.getDerivativeIncrementLength() * normx
643 else:
644 epsnew = self.getDerivativeIncrementLength()
645 if not v_is_normalised:
646 normv=self.norm(v)
647 if normv<=0:
648 return F0*0
649 else:
650 epsnew /= normv
651 F1=self.eval(x0 + epsnew * v)
652 return (F1-F0)/epsnew
653
654 ######################################
655 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):
656 """
657 solves a non-linear problem M{F(x)=0} for unknown M{x} using the stopping criterion:
658
659 M{norm(F(x) <= atol + rtol * norm(F(x0)}
660
661 where M{x0} is the initial guess.
662
663 @param defect: object defining the the function M{F}, C{defect.norm} defines the M{norm} used in the stopping criterion.
664 @type defect: L{Defect}
665 @param x: initial guess for the solution, C{x} is altered.
666 @type x: any object type allowing basic operations such as L{numarray.NumArray}, L{Data}
667 @param iter_max: maximum number of iteration steps
668 @type iter_max: positive C{int}
669 @param sub_iter_max:
670 @type sub_iter_max:
671 @param atol: absolute tolerance for the solution
672 @type atol: positive C{float}
673 @param rtol: relative tolerance for the solution
674 @type rtol: positive C{float}
675 @param gamma: tolerance safety factor for inner iteration
676 @type gamma: positive C{float}, less than 1
677 @param sub_tol_max: upper bound for inner tolerance.
678 @type sub_tol_max: positive C{float}, less than 1
679 @return: an approximation of the solution with the desired accuracy
680 @rtype: same type as the initial guess.
681 """
682 lmaxit=iter_max
683 if atol<0: raise ValueError,"atol needs to be non-negative."
684 if rtol<0: raise ValueError,"rtol needs to be non-negative."
685 if rtol+atol<=0: raise ValueError,"rtol or atol needs to be non-negative."
686 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
687 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
688
689 F=defect(x)
690 fnrm=defect.norm(F)
691 stop_tol=atol + rtol*fnrm
692 sub_tol=sub_tol_max
693 if verbose: print "NewtonGMRES: initial residual = %e."%fnrm
694 if verbose: print " tolerance = %e."%sub_tol
695 iter=1
696 #
697 # main iteration loop
698 #
699 while not fnrm<=stop_tol:
700 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
701 #
702 # adjust sub_tol_
703 #
704 if iter > 1:
705 rat=fnrm/fnrmo
706 sub_tol_old=sub_tol
707 sub_tol=gamma*rat**2
708 if gamma*sub_tol_old**2 > .1: sub_tol=max(sub_tol,gamma*sub_tol_old**2)
709 sub_tol=max(min(sub_tol,sub_tol_max), .5*stop_tol/fnrm)
710 #
711 # calculate newton increment xc
712 # if iter_max in __FDGMRES is reached MaxIterReached is thrown
713 # if iter_restart -1 is returned as sub_iter
714 # if atol is reached sub_iter returns the numer of steps performed to get there
715 #
716 #
717 if verbose: print " subiteration (GMRES) is called with relative tolerance %e."%sub_tol
718 try:
719 xc, sub_iter=__FDGMRES(F, defect, x, sub_tol*fnrm, iter_max=iter_max-iter, iter_restart=sub_iter_max)
720 except MaxIterReached:
721 raise MaxIterReached,"maximum number of %s steps reached."%iter_max
722 if sub_iter<0:
723 iter+=sub_iter_max
724 else:
725 iter+=sub_iter
726 # ====
727 x+=xc
728 F=defect(x)
729 iter+=1
730 fnrmo, fnrm=fnrm, defect.norm(F)
731 if verbose: print " step %s: residual %e."%(iter,fnrm)
732 if verbose: print "NewtonGMRES: completed after %s steps."%iter
733 return x
734
735 def __givapp(c,s,vin):
736 """
737 apply a sequence of Givens rotations (c,s) to the recuirsively to the vector vin
738 @warning: C{vin} is altered.
739 """
740 vrot=vin
741 if isinstance(c,float):
742 vrot=[c*vrot[0]-s*vrot[1],s*vrot[0]+c*vrot[1]]
743 else:
744 for i in range(len(c)):
745 w1=c[i]*vrot[i]-s[i]*vrot[i+1]
746 w2=s[i]*vrot[i]+c[i]*vrot[i+1]
747 vrot[i:i+2]=w1,w2
748 return vrot
749
750 def __FDGMRES(F0, defect, x0, atol, iter_max=100, iter_restart=20):
751 h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)
752 c=numarray.zeros(iter_restart,numarray.Float64)
753 s=numarray.zeros(iter_restart,numarray.Float64)
754 g=numarray.zeros(iter_restart,numarray.Float64)
755 v=[]
756
757 rho=defect.norm(F0)
758 if rho<=0.: return x0*0
759
760 v.append(-F0/rho)
761 g[0]=rho
762 iter=0
763 while rho > atol and iter<iter_restart-1:
764
765 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
766
767 p=defect.derivative(F0,x0,v[iter], v_is_normalised=True)
768 v.append(p)
769
770 v_norm1=defect.norm(v[iter+1])
771
772 # Modified Gram-Schmidt
773 for j in range(iter+1):
774 h[j][iter]=defect.bilinearform(v[j],v[iter+1])
775 v[iter+1]-=h[j][iter]*v[j]
776
777 h[iter+1][iter]=defect.norm(v[iter+1])
778 v_norm2=h[iter+1][iter]
779
780 # Reorthogonalize if needed
781 if v_norm1 + 0.001*v_norm2 == v_norm1: #Brown/Hindmarsh condition (default)
782 for j in range(iter+1):
783 hr=defect.bilinearform(v[j],v[iter+1])
784 h[j][iter]=h[j][iter]+hr
785 v[iter+1] -= hr*v[j]
786
787 v_norm2=defect.norm(v[iter+1])
788 h[iter+1][iter]=v_norm2
789 # watch out for happy breakdown
790 if not v_norm2 == 0:
791 v[iter+1]=v[iter+1]/h[iter+1][iter]
792
793 # Form and store the information for the new Givens rotation
794 if iter > 0 :
795 hhat=numarray.zeros(iter+1,numarray.Float64)
796 for i in range(iter+1) : hhat[i]=h[i][iter]
797 hhat=__givapp(c[0:iter],s[0:iter],hhat);
798 for i in range(iter+1) : h[i][iter]=hhat[i]
799
800 mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])
801
802 if mu!=0 :
803 c[iter]=h[iter][iter]/mu
804 s[iter]=-h[iter+1][iter]/mu
805 h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter]
806 h[iter+1][iter]=0.0
807 g[iter:iter+2]=__givapp(c[iter],s[iter],g[iter:iter+2])
808
809 # Update the residual norm
810 rho=abs(g[iter+1])
811 iter+=1
812
813 # At this point either iter > iter_max or rho < tol.
814 # It's time to compute x and leave.
815 if iter > 0 :
816 y=numarray.zeros(iter,numarray.Float64)
817 y[iter-1] = g[iter-1] / h[iter-1][iter-1]
818 if iter > 1 :
819 i=iter-2
820 while i>=0 :
821 y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]
822 i=i-1
823 xhat=v[iter-1]*y[iter-1]
824 for i in range(iter-1):
825 xhat += v[i]*y[i]
826 else :
827 xhat=v[0] * 0
828
829 if iter<iter_restart-1:
830 stopped=iter
831 else:
832 stopped=-1
833
834 return xhat,stopped
835
836 ##############################################
837 def GMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):
838 ################################################
839 m=iter_restart
840 iter=0
841 xc=x
842 while True:
843 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached"%iter_max
844 xc,stopped=__GMRESm(b*1, Aprod, Msolve, bilinearform, stoppingcriterium, x=xc*1, iter_max=iter_max-iter, iter_restart=m)
845 if stopped: break
846 iter+=iter_restart
847 return xc
848
849 def __GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):
850 iter=0
851 r=Msolve(b)
852 r_dot_r = bilinearform(r, r)
853 if r_dot_r<0: raise NegativeNorm,"negative norm."
854 norm_b=math.sqrt(r_dot_r)
855
856 if x==None:
857 x=0*b
858 else:
859 r=Msolve(b-Aprod(x))
860 r_dot_r = bilinearform(r, r)
861 if r_dot_r<0: raise NegativeNorm,"negative norm."
862
863 h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)
864 c=numarray.zeros(iter_restart,numarray.Float64)
865 s=numarray.zeros(iter_restart,numarray.Float64)
866 g=numarray.zeros(iter_restart,numarray.Float64)
867 v=[]
868
869 rho=math.sqrt(r_dot_r)
870
871 v.append(r/rho)
872 g[0]=rho
873
874 while not (stoppingcriterium(rho,norm_b) or iter==iter_restart-1):
875
876 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
877
878 p=Msolve(Aprod(v[iter]))
879
880 v.append(p)
881
882 v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
883
884 # Modified Gram-Schmidt
885 for j in range(iter+1):
886 h[j][iter]=bilinearform(v[j],v[iter+1])
887 v[iter+1]-=h[j][iter]*v[j]
888
889 h[iter+1][iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
890 v_norm2=h[iter+1][iter]
891
892 # Reorthogonalize if needed
893 if v_norm1 + 0.001*v_norm2 == v_norm1: #Brown/Hindmarsh condition (default)
894 for j in range(iter+1):
895 hr=bilinearform(v[j],v[iter+1])
896 h[j][iter]=h[j][iter]+hr
897 v[iter+1] -= hr*v[j]
898
899 v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
900 h[iter+1][iter]=v_norm2
901
902 # watch out for happy breakdown
903 if not v_norm2 == 0:
904 v[iter+1]=v[iter+1]/h[iter+1][iter]
905
906 # Form and store the information for the new Givens rotation
907 if iter > 0 :
908 hhat=numarray.zeros(iter+1,numarray.Float64)
909 for i in range(iter+1) : hhat[i]=h[i][iter]
910 hhat=__givapp(c[0:iter],s[0:iter],hhat);
911 for i in range(iter+1) : h[i][iter]=hhat[i]
912
913 mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])
914
915 if mu!=0 :
916 c[iter]=h[iter][iter]/mu
917 s[iter]=-h[iter+1][iter]/mu
918 h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter]
919 h[iter+1][iter]=0.0
920 g[iter:iter+2]=__givapp(c[iter],s[iter],g[iter:iter+2])
921
922 # Update the residual norm
923
924 rho=abs(g[iter+1])
925 iter+=1
926
927 # At this point either iter > iter_max or rho < tol.
928 # It's time to compute x and leave.
929
930 if iter > 0 :
931 y=numarray.zeros(iter,numarray.Float64)
932 y[iter-1] = g[iter-1] / h[iter-1][iter-1]
933 if iter > 1 :
934 i=iter-2
935 while i>=0 :
936 y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]
937 i=i-1
938 xhat=v[iter-1]*y[iter-1]
939 for i in range(iter-1):
940 xhat += v[i]*y[i]
941 else : xhat=v[0]
942
943 x += xhat
944 if iter<iter_restart-1:
945 stopped=True
946 else:
947 stopped=False
948
949 return x,stopped
950
951 #################################################
952 def MINRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
953 #################################################
954 #
955 # minres solves the system of linear equations Ax = b
956 # where A is a symmetric matrix (possibly indefinite or singular)
957 # and b is a given vector.
958 #
959 # "A" may be a dense or sparse matrix (preferably sparse!)
960 # or the name of a function such that
961 # y = A(x)
962 # returns the product y = Ax for any given vector x.
963 #
964 # "M" defines a positive-definite preconditioner M = C C'.
965 # "M" may be a dense or sparse matrix (preferably sparse!)
966 # or the name of a function such that
967 # solves the system My = x for any given vector x.
968 #
969 #
970
971 #------------------------------------------------------------------
972 # Set up y and v for the first Lanczos vector v1.
973 # y = beta1 P' v1, where P = C**(-1).
974 # v is really P' v1.
975 #------------------------------------------------------------------
976 if x==None:
977 x=0*b
978 else:
979 b += (-1)*Aprod(x)
980
981 r1 = b
982 y = Msolve(b)
983 beta1 = bilinearform(y,b)
984
985 if beta1< 0: raise NegativeNorm,"negative norm."
986
987 # If b = 0 exactly, stop with x = 0.
988 if beta1==0: return x*0.
989
990 if beta1> 0:
991 beta1 = math.sqrt(beta1)
992
993 #------------------------------------------------------------------
994 # Initialize quantities.
995 # ------------------------------------------------------------------
996 iter = 0
997 Anorm = 0
998 ynorm = 0
999 oldb = 0
1000 beta = beta1
1001 dbar = 0
1002 epsln = 0
1003 phibar = beta1
1004 rhs1 = beta1
1005 rhs2 = 0
1006 rnorm = phibar
1007 tnorm2 = 0
1008 ynorm2 = 0
1009 cs = -1
1010 sn = 0
1011 w = b*0.
1012 w2 = b*0.
1013 r2 = r1
1014 eps = 0.0001
1015
1016 #---------------------------------------------------------------------
1017 # Main iteration loop.
1018 # --------------------------------------------------------------------
1019 while not stoppingcriterium(rnorm,Anorm*ynorm,'MINRES'): # checks ||r|| < (||A|| ||x||) * TOL
1020
1021 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
1022 iter = iter + 1
1023
1024 #-----------------------------------------------------------------
1025 # Obtain quantities for the next Lanczos vector vk+1, k = 1, 2,...
1026 # The general iteration is similar to the case k = 1 with v0 = 0:
1027 #
1028 # p1 = Operator * v1 - beta1 * v0,
1029 # alpha1 = v1'p1,
1030 # q2 = p2 - alpha1 * v1,
1031 # beta2^2 = q2'q2,
1032 # v2 = (1/beta2) q2.
1033 #
1034 # Again, y = betak P vk, where P = C**(-1).
1035 #-----------------------------------------------------------------
1036 s = 1/beta # Normalize previous vector (in y).
1037 v = s*y # v = vk if P = I
1038
1039 y = Aprod(v)
1040
1041 if iter >= 2:
1042 y = y - (beta/oldb)*r1
1043
1044 alfa = bilinearform(v,y) # alphak
1045 y += (- alfa/beta)*r2
1046 r1 = r2
1047 r2 = y
1048 y = Msolve(r2)
1049 oldb = beta # oldb = betak
1050 beta = bilinearform(y,r2) # beta = betak+1^2
1051 if beta < 0: raise NegativeNorm,"negative norm."
1052
1053 beta = math.sqrt( beta )
1054 tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta
1055
1056 if iter==1: # Initialize a few things.
1057 gmax = abs( alfa ) # alpha1
1058 gmin = gmax # alpha1
1059
1060 # Apply previous rotation Qk-1 to get
1061 # [deltak epslnk+1] = [cs sn][dbark 0 ]
1062 # [gbar k dbar k+1] [sn -cs][alfak betak+1].
1063
1064 oldeps = epsln
1065 delta = cs * dbar + sn * alfa # delta1 = 0 deltak
1066 gbar = sn * dbar - cs * alfa # gbar 1 = alfa1 gbar k
1067 epsln = sn * beta # epsln2 = 0 epslnk+1
1068 dbar = - cs * beta # dbar 2 = beta2 dbar k+1
1069
1070 # Compute the next plane rotation Qk
1071
1072 gamma = math.sqrt(gbar*gbar+beta*beta) # gammak
1073 gamma = max(gamma,eps)
1074 cs = gbar / gamma # ck
1075 sn = beta / gamma # sk
1076 phi = cs * phibar # phik
1077 phibar = sn * phibar # phibark+1
1078
1079 # Update x.
1080
1081 denom = 1/gamma
1082 w1 = w2
1083 w2 = w
1084 w = (v - oldeps*w1 - delta*w2) * denom
1085 x += phi*w
1086
1087 # Go round again.
1088
1089 gmax = max(gmax,gamma)
1090 gmin = min(gmin,gamma)
1091 z = rhs1 / gamma
1092 ynorm2 = z*z + ynorm2
1093 rhs1 = rhs2 - delta*z
1094 rhs2 = - epsln*z
1095
1096 # Estimate various norms and test for convergence.
1097
1098 Anorm = math.sqrt( tnorm2 )
1099 ynorm = math.sqrt( ynorm2 )
1100
1101 rnorm = phibar
1102
1103 return x
1104
1105 def TFQMR(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
1106
1107 # TFQMR solver for linear systems
1108 #
1109 #
1110 # initialization
1111 #
1112 errtol = math.sqrt(bilinearform(b,b))
1113 norm_b=errtol
1114 kmax = iter_max
1115 error = []
1116
1117 if math.sqrt(bilinearform(x,x)) != 0.0:
1118 r = b - Aprod(x)
1119 else:
1120 r = b
1121
1122 r=Msolve(r)
1123
1124 u1=0
1125 u2=0
1126 y1=0
1127 y2=0
1128
1129 w = r
1130 y1 = r
1131 iter = 0
1132 d = 0
1133
1134 v = Msolve(Aprod(y1))
1135 u1 = v
1136
1137 theta = 0.0;
1138 eta = 0.0;
1139 tau = math.sqrt(bilinearform(r,r))
1140 error = [ error, tau ]
1141 rho = tau * tau
1142 m=1
1143 #
1144 # TFQMR iteration
1145 #
1146 # while ( iter < kmax-1 ):
1147
1148 while not stoppingcriterium(tau*math.sqrt ( m + 1 ),norm_b,'TFQMR'):
1149 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
1150
1151 sigma = bilinearform(r,v)
1152
1153 if ( sigma == 0.0 ):
1154 raise 'TFQMR breakdown, sigma=0'
1155
1156
1157 alpha = rho / sigma
1158
1159 for j in range(2):
1160 #
1161 # Compute y2 and u2 only if you have to
1162 #
1163 if ( j == 1 ):
1164 y2 = y1 - alpha * v
1165 u2 = Msolve(Aprod(y2))
1166
1167 m = 2 * (iter+1) - 2 + (j+1)
1168 if j==0:
1169 w = w - alpha * u1
1170 d = y1 + ( theta * theta * eta / alpha ) * d
1171 if j==1:
1172 w = w - alpha * u2
1173 d = y2 + ( theta * theta * eta / alpha ) * d
1174
1175 theta = math.sqrt(bilinearform(w,w))/ tau
1176 c = 1.0 / math.sqrt ( 1.0 + theta * theta )
1177 tau = tau * theta * c
1178 eta = c * c * alpha
1179 x = x + eta * d
1180 #
1181 # Try to terminate the iteration at each pass through the loop
1182 #
1183 # if ( tau * math.sqrt ( m + 1 ) <= errtol ):
1184 # error = [ error, tau ]
1185 # total_iters = iter
1186 # break
1187
1188
1189 if ( rho == 0.0 ):
1190 raise 'TFQMR breakdown, rho=0'
1191
1192
1193 rhon = bilinearform(r,w)
1194 beta = rhon / rho;
1195 rho = rhon;
1196 y1 = w + beta * y2;
1197 u1 = Msolve(Aprod(y1))
1198 v = u1 + beta * ( u2 + beta * v )
1199 error = [ error, tau ]
1200 total_iters = iter
1201
1202 iter = iter + 1
1203
1204 return x
1205
1206
1207 #############################################
1208
1209 class ArithmeticTuple(object):
1210 """
1211 tuple supporting inplace update x+=y and scaling x=a*y where x,y is an ArithmeticTuple and a is a float.
1212
1213 example of usage:
1214
1215 from esys.escript import Data
1216 from numarray import array
1217 a=Data(...)
1218 b=array([1.,4.])
1219 x=ArithmeticTuple(a,b)
1220 y=5.*x
1221
1222 """
1223 def __init__(self,*args):
1224 """
1225 initialize object with elements args.
1226
1227 @param args: tuple of object that support implace add (x+=y) and scaling (x=a*y)
1228 """
1229 self.__items=list(args)
1230
1231 def __len__(self):
1232 """
1233 number of items
1234
1235 @return: number of items
1236 @rtype: C{int}
1237 """
1238 return len(self.__items)
1239
1240 def __getitem__(self,index):
1241 """
1242 get an item
1243
1244 @param index: item to be returned
1245 @type index: C{int}
1246 @return: item with index C{index}
1247 """
1248 return self.__items.__getitem__(index)
1249
1250 def __mul__(self,other):
1251 """
1252 scaling from the right
1253
1254 @param other: scaling factor
1255 @type other: C{float}
1256 @return: itemwise self*other
1257 @rtype: L{ArithmeticTuple}
1258 """
1259 out=[]
1260 try:
1261 l=len(other)
1262 if l!=len(self):
1263 raise ValueError,"length of of arguments don't match."
1264 for i in range(l): out.append(self[i]*other[i])
1265 except TypeError:
1266 for i in range(len(self)): out.append(self[i]*other)
1267 return ArithmeticTuple(*tuple(out))
1268
1269 def __rmul__(self,other):
1270 """
1271 scaling from the left
1272
1273 @param other: scaling factor
1274 @type other: C{float}
1275 @return: itemwise other*self
1276 @rtype: L{ArithmeticTuple}
1277 """
1278 out=[]
1279 try:
1280 l=len(other)
1281 if l!=len(self):
1282 raise ValueError,"length of of arguments don't match."
1283 for i in range(l): out.append(other[i]*self[i])
1284 except TypeError:
1285 for i in range(len(self)): out.append(other*self[i])
1286 return ArithmeticTuple(*tuple(out))
1287
1288 def __div__(self,other):
1289 """
1290 dividing from the right
1291
1292 @param other: scaling factor
1293 @type other: C{float}
1294 @return: itemwise self/other
1295 @rtype: L{ArithmeticTuple}
1296 """
1297 return self*(1/other)
1298
1299 def __rdiv__(self,other):
1300 """
1301 dividing from the left
1302
1303 @param other: scaling factor
1304 @type other: C{float}
1305 @return: itemwise other/self
1306 @rtype: L{ArithmeticTuple}
1307 """
1308 out=[]
1309 try:
1310 l=len(other)
1311 if l!=len(self):
1312 raise ValueError,"length of of arguments don't match."
1313 for i in range(l): out.append(other[i]/self[i])
1314 except TypeError:
1315 for i in range(len(self)): out.append(other/self[i])
1316 return ArithmeticTuple(*tuple(out))
1317
1318 def __iadd__(self,other):
1319 """
1320 in-place add of other to self
1321
1322 @param other: increment
1323 @type other: C{ArithmeticTuple}
1324 """
1325 if len(self) != len(other):
1326 raise ValueError,"tuple length must match."
1327 for i in range(len(self)):
1328 self.__items[i]+=other[i]
1329 return self
1330
1331 def __add__(self,other):
1332 """
1333 add other to self
1334
1335 @param other: increment
1336 @type other: C{ArithmeticTuple}
1337 """
1338 out=[]
1339 try:
1340 l=len(other)
1341 if l!=len(self):
1342 raise ValueError,"length of of arguments don't match."
1343 for i in range(l): out.append(self[i]+other[i])
1344 except TypeError:
1345 for i in range(len(self)): out.append(self[i]+other)
1346 return ArithmeticTuple(*tuple(out))
1347
1348 def __sub__(self,other):
1349 """
1350 subtract other from self
1351
1352 @param other: increment
1353 @type other: C{ArithmeticTuple}
1354 """
1355 out=[]
1356 try:
1357 l=len(other)
1358 if l!=len(self):
1359 raise ValueError,"length of of arguments don't match."
1360 for i in range(l): out.append(self[i]-other[i])
1361 except TypeError:
1362 for i in range(len(self)): out.append(self[i]-other)
1363 return ArithmeticTuple(*tuple(out))
1364
1365 def __isub__(self,other):
1366 """
1367 subtract other from self
1368
1369 @param other: increment
1370 @type other: C{ArithmeticTuple}
1371 """
1372 if len(self) != len(other):
1373 raise ValueError,"tuple length must match."
1374 for i in range(len(self)):
1375 self.__items[i]-=other[i]
1376 return self
1377
1378 def __neg__(self):
1379 """
1380 negate
1381
1382 """
1383 out=[]
1384 for i in range(len(self)):
1385 out.append(-self[i])
1386 return ArithmeticTuple(*tuple(out))
1387
1388
1389 class HomogeneousSaddlePointProblem(object):
1390 """
1391 This provides a framwork for solving linear homogeneous saddle point problem of the form
1392
1393 Av+B^*p=f
1394 Bv =0
1395
1396 for the unknowns v and p and given operators A and B and given right hand side f.
1397 B^* is the adjoint operator of B is the given inner product.
1398
1399 """
1400 def __init__(self,**kwargs):
1401 self.setTolerance()
1402 self.setToleranceReductionFactor()
1403
1404 def initialize(self):
1405 """
1406 initialize the problem (overwrite)
1407 """
1408 pass
1409 def B(self,v):
1410 """
1411 returns Bv (overwrite)
1412 @rtype: equal to the type of p
1413
1414 @note: boundary conditions on p should be zero!
1415 """
1416 pass
1417
1418 def inner(self,p0,p1):
1419 """
1420 returns inner product of two element p0 and p1 (overwrite)
1421
1422 @type p0: equal to the type of p
1423 @type p1: equal to the type of p
1424 @rtype: C{float}
1425
1426 @rtype: equal to the type of p
1427 """
1428 pass
1429
1430 def solve_A(self,u,p):
1431 """
1432 solves Av=f-Au-B^*p with accuracy self.getReducedTolerance() (overwrite)
1433
1434 @rtype: equal to the type of v
1435 @note: boundary conditions on v should be zero!
1436 """
1437 pass
1438
1439 def solve_prec(self,p):
1440 """
1441 provides a preconditioner for BA^{-1}B^* with accuracy self.getReducedTolerance() (overwrite)
1442
1443 @rtype: equal to the type of p
1444 """
1445 pass
1446
1447 def stoppingcriterium(self,Bv,v,p):
1448 """
1449 returns a True if iteration is terminated. (overwrite)
1450
1451 @rtype: C{bool}
1452 """
1453 pass
1454
1455 def __inner(self,p,r):
1456 return self.inner(p,r[1])
1457
1458 def __inner_p(self,p1,p2):
1459 return self.inner(p1,p2)
1460
1461 def __inner_a(self,a1,a2):
1462 return self.inner_a(a1,a2)
1463
1464 def __inner_a1(self,a1,a2):
1465 return self.inner(a1[1],a2[1])
1466
1467 def __stoppingcriterium(self,norm_r,r,p):
1468 return self.stoppingcriterium(r[1],r[0],p)
1469
1470 def __stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None):
1471 return self.stoppingcriterium2(norm_r,norm_b,solver,TOL)
1472
1473 def setTolerance(self,tolerance=1.e-8):
1474 self.__tol=tolerance
1475 def getTolerance(self):
1476 return self.__tol
1477 def setToleranceReductionFactor(self,reduction=0.01):
1478 self.__reduction=reduction
1479 def getSubProblemTolerance(self):
1480 return self.__reduction*self.getTolerance()
1481
1482 def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='PCG',iter_restart=20):
1483 """
1484 solves the saddle point problem using initial guesses v and p.
1485
1486 @param max_iter: maximum number of iteration steps.
1487 """
1488 self.verbose=verbose
1489 self.show_details=show_details and self.verbose
1490
1491 # assume p is known: then v=A^-1(f-B^*p)
1492 # which leads to BA^-1B^*p = BA^-1f
1493
1494 # Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask)
1495 self.__z=v+self.solve_A(v,p*0)
1496 Bz=self.B(self.__z)
1497 #
1498 # solve BA^-1B^*p = Bz
1499 #
1500 #
1501 #
1502 self.iter=0
1503 if solver=='GMRES':
1504 if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter
1505 p=GMRES(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.,iter_restart=iter_restart)
1506 # solve Au=f-B^*p
1507 # A(u-v)=f-B^*p-Av
1508 # u=v+(u-v)
1509 u=v+self.solve_A(v,p)
1510
1511 if solver=='TFQMR':
1512 if self.verbose: print "enter TFQMR method (iter_max=%s)"%max_iter
1513 p=TFQMR(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)
1514 # solve Au=f-B^*p
1515 # A(u-v)=f-B^*p-Av
1516 # u=v+(u-v)
1517 u=v+self.solve_A(v,p)
1518
1519 if solver=='MINRES':
1520 if self.verbose: print "enter MINRES method (iter_max=%s)"%max_iter
1521 p=MINRES(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)
1522 # solve Au=f-B^*p
1523 # A(u-v)=f-B^*p-Av
1524 # u=v+(u-v)
1525 u=v+self.solve_A(v,p)
1526
1527 if solver=='GMRESC':
1528 if self.verbose: print "enter GMRES coupled method (iter_max=%s)"%max_iter
1529 p0=self.solve_prec1(Bz)
1530 #alfa=(1/self.vol)*util.integrate(util.interpolate(p,escript.Function(self.domain)))
1531 #p-=alfa
1532 x=GMRES(ArithmeticTuple(self.__z*1.,p0*1),self.__Anext,self.__Mempty,self.__inner_a,self.__stoppingcriterium2,iter_max=max_iter, x=ArithmeticTuple(v*1,p*1),iter_restart=20)
1533 #x=NewtonGMRES(ArithmeticTuple(self.__z*1.,p0*1),self.__Aprod_Newton2,self.__Mempty,self.__inner_a,self.__stoppingcriterium2,iter_max=max_iter, x=ArithmeticTuple(v*1,p*1),atol=0,rtol=self.getTolerance())
1534
1535 # solve Au=f-B^*p
1536 # A(u-v)=f-B^*p-Av
1537 # u=v+(u-v)
1538 p=x[1]
1539 u=v+self.solve_A(v,p)
1540 #p=x[1]
1541 #u=x[0]
1542
1543 if solver=='PCG':
1544 # note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv
1545 #
1546 # with Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask)
1547 # A(v-z)= f -Az - B^*p (v-z=0 on fixed_u_mask)
1548 if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter
1549 p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p)
1550 u=r[0]
1551 # print "DIFF=",util.integrate(p)
1552
1553 # print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)
1554
1555 return u,p
1556
1557 def __Msolve(self,r):
1558 return self.solve_prec(r[1])
1559
1560 def __Msolve2(self,r):
1561 return self.solve_prec(r*1)
1562
1563 def __Mempty(self,r):
1564 return r
1565
1566
1567 def __Aprod(self,p):
1568 # return BA^-1B*p
1569 #solve Av =B^*p as Av =f-Az-B^*(-p)
1570 v=self.solve_A(self.__z,-p)
1571 return ArithmeticTuple(v, self.B(v))
1572
1573 def __Anext(self,x):
1574 # return next v,p
1575 #solve Av =-B^*p as Av =f-Az-B^*p
1576
1577 pc=x[1]
1578 v=self.solve_A(self.__z,-pc)
1579 p=self.solve_prec1(self.B(v))
1580
1581 return ArithmeticTuple(v,p)
1582
1583
1584 def __Aprod2(self,p):
1585 # return BA^-1B*p
1586 #solve Av =B^*p as Av =f-Az-B^*(-p)
1587 v=self.solve_A(self.__z,-p)
1588 return self.B(v)
1589
1590 def __Aprod_Newton(self,p):
1591 # return BA^-1B*p - Bz
1592 #solve Av =-B^*p as Av =f-Az-B^*p
1593 v=self.solve_A(self.__z,-p)
1594 return self.B(v-self.__z)
1595
1596 def __Aprod_Newton2(self,x):
1597 # return BA^-1B*p - Bz
1598 #solve Av =-B^*p as Av =f-Az-B^*p
1599 pc=x[1]
1600 v=self.solve_A(self.__z,-pc)
1601 p=self.solve_prec1(self.B(v-self.__z))
1602 return ArithmeticTuple(v,p)
1603
1604
1605 def MaskFromBoundaryTag(domain,*tags):
1606 """
1607 create a mask on the Solution(domain) function space which one for samples
1608 that touch the boundary tagged by tags.
1609
1610 usage: m=MaskFromBoundaryTag(domain,"left", "right")
1611
1612 @param domain: a given domain
1613 @type domain: L{escript.Domain}
1614 @param tags: boundray tags
1615 @type tags: C{str}
1616 @return: a mask which marks samples that are touching the boundary tagged by any of the given tags.
1617 @rtype: L{escript.Data} of rank 0
1618 """
1619 pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1620 d=escript.Scalar(0.,escript.FunctionOnBoundary(function_space.getDomain()))
1621 for t in tags: d.setTaggedValue(t,1.)
1622 pde.setValue(y=d)
1623 return util.whereNonZero(pde.getRightHandSide())
1624 #============================================================================================================================================
1625 class SaddlePointProblem(object):
1626 """
1627 This implements a solver for a saddlepoint problem
1628
1629 M{f(u,p)=0}
1630 M{g(u)=0}
1631
1632 for u and p. The problem is solved with an inexact Uszawa scheme for p:
1633
1634 M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})}
1635 M{Q_g (p^{k+1}-p^{k}) = g(u^{k+1})}
1636
1637 where Q_f is an approximation of the Jacobiean A_f of f with respect to u and Q_f is an approximation of
1638 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'
1639 Q_g can be difficult, non-linear conjugate gradient method is applied to solve for p, so Q_g plays
1640 in fact the role of a preconditioner.
1641 """
1642 def __init__(self,verbose=False,*args):
1643 """
1644 initializes the problem
1645
1646 @param verbose: switches on the printing out some information
1647 @type verbose: C{bool}
1648 @note: this method may be overwritten by a particular saddle point problem
1649 """
1650 print "SaddlePointProblem should not be used anymore!"
1651 if not isinstance(verbose,bool):
1652 raise TypeError("verbose needs to be of type bool.")
1653 self.__verbose=verbose
1654 self.relaxation=1.
1655 DeprecationWarning("SaddlePointProblem should not be used anymore.")
1656
1657 def trace(self,text):
1658 """
1659 prints text if verbose has been set
1660
1661 @param text: a text message
1662 @type text: C{str}
1663 """
1664 if self.__verbose: print "%s: %s"%(str(self),text)
1665
1666 def solve_f(self,u,p,tol=1.e-8):
1667 """
1668 solves
1669
1670 A_f du = f(u,p)
1671
1672 with tolerance C{tol} and return du. A_f is Jacobiean of f with respect to u.
1673
1674 @param u: current approximation of u
1675 @type u: L{escript.Data}
1676 @param p: current approximation of p
1677 @type p: L{escript.Data}
1678 @param tol: tolerance expected for du
1679 @type tol: C{float}
1680 @return: increment du
1681 @rtype: L{escript.Data}
1682 @note: this method has to be overwritten by a particular saddle point problem
1683 """
1684 pass
1685
1686 def solve_g(self,u,tol=1.e-8):
1687 """
1688 solves
1689
1690 Q_g dp = g(u)
1691
1692 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.
1693
1694 @param u: current approximation of u
1695 @type u: L{escript.Data}
1696 @param tol: tolerance expected for dp
1697 @type tol: C{float}
1698 @return: increment dp
1699 @rtype: L{escript.Data}
1700 @note: this method has to be overwritten by a particular saddle point problem
1701 """
1702 pass
1703
1704 def inner(self,p0,p1):
1705 """
1706 inner product of p0 and p1 approximating p. Typically this returns integrate(p0*p1)
1707 @return: inner product of p0 and p1
1708 @rtype: C{float}
1709 """
1710 pass
1711
1712 subiter_max=3
1713 def solve(self,u0,p0,tolerance=1.e-6,tolerance_u=None,iter_max=100,accepted_reduction=0.995,relaxation=None):
1714 """
1715 runs the solver
1716
1717 @param u0: initial guess for C{u}
1718 @type u0: L{esys.escript.Data}
1719 @param p0: initial guess for C{p}
1720 @type p0: L{esys.escript.Data}
1721 @param tolerance: tolerance for relative error in C{u} and C{p}
1722 @type tolerance: positive C{float}
1723 @param tolerance_u: tolerance for relative error in C{u} if different from C{tolerance}
1724 @type tolerance_u: positive C{float}
1725 @param iter_max: maximum number of iteration steps.
1726 @type iter_max: C{int}
1727 @param accepted_reduction: if the norm g cannot be reduced by C{accepted_reduction} backtracking to adapt the
1728 relaxation factor. If C{accepted_reduction=None} no backtracking is used.
1729 @type accepted_reduction: positive C{float} or C{None}
1730 @param relaxation: initial relaxation factor. If C{relaxation==None}, the last relaxation factor is used.
1731 @type relaxation: C{float} or C{None}
1732 """
1733 tol=1.e-2
1734 if tolerance_u==None: tolerance_u=tolerance
1735 if not relaxation==None: self.relaxation=relaxation
1736 if accepted_reduction ==None:
1737 angle_limit=0.
1738 elif accepted_reduction>=1.:
1739 angle_limit=0.
1740 else:
1741 angle_limit=util.sqrt(1-accepted_reduction**2)
1742 self.iter=0
1743 u=u0
1744 p=p0
1745 #
1746 # initialize things:
1747 #
1748 converged=False
1749 #
1750 # start loop:
1751 #
1752 # initial search direction is g
1753 #
1754 while not converged :
1755 if self.iter>iter_max:
1756 raise ArithmeticError("no convergence after %s steps."%self.iter)
1757 f_new=self.solve_f(u,p,tol)
1758 norm_f_new = util.Lsup(f_new)
1759 u_new=u-f_new
1760 g_new=self.solve_g(u_new,tol)
1761 self.iter+=1
1762 norm_g_new = util.sqrt(self.inner(g_new,g_new))
1763 if norm_f_new==0. and norm_g_new==0.: return u, p
1764 if self.iter>1 and not accepted_reduction==None:
1765 #
1766 # did we manage to reduce the norm of G? I
1767 # if not we start a backtracking procedure
1768 #
1769 # print "new/old norm = ",norm_g_new, norm_g, norm_g_new/norm_g
1770 if norm_g_new > accepted_reduction * norm_g:
1771 sub_iter=0
1772 s=self.relaxation
1773 d=g
1774 g_last=g
1775 self.trace(" start substepping: f = %s, g = %s, relaxation = %s."%(norm_f_new, norm_g_new, s))
1776 while sub_iter < self.subiter_max and norm_g_new > accepted_reduction * norm_g:
1777 dg= g_new-g_last
1778 norm_dg=abs(util.sqrt(self.inner(dg,dg))/self.relaxation)
1779 rad=self.inner(g_new,dg)/self.relaxation
1780 # print " ",sub_iter,": rad, norm_dg:",abs(rad), norm_dg*norm_g_new * angle_limit
1781 # print " ",sub_iter,": rad, norm_dg:",rad, norm_dg, norm_g_new, norm_g
1782 if abs(rad) < norm_dg*norm_g_new * angle_limit:
1783 if sub_iter>0: self.trace(" no further improvements expected from backtracking.")
1784 break
1785 r=self.relaxation
1786 self.relaxation= - rad/norm_dg**2
1787 s+=self.relaxation
1788 #####
1789 # a=g_new+self.relaxation*dg/r
1790 # print "predicted new norm = ",util.sqrt(self.inner(a,a)),util.sqrt(self.inner(g_new,g_new)), self.relaxation
1791 #####
1792 g_last=g_new
1793 p+=self.relaxation*d
1794 f_new=self.solve_f(u,p,tol)
1795 u_new=u-f_new
1796 g_new=self.solve_g(u_new,tol)
1797 self.iter+=1
1798 norm_f_new = util.Lsup(f_new)
1799 norm_g_new = util.sqrt(self.inner(g_new,g_new))
1800 # print " ",sub_iter," new g norm",norm_g_new
1801 self.trace(" %s th sub-step: f = %s, g = %s, relaxation = %s."%(sub_iter, norm_f_new, norm_g_new, s))
1802 #
1803 # can we expect reduction of g?
1804 #
1805 # u_last=u_new
1806 sub_iter+=1
1807 self.relaxation=s
1808 #
1809 # check for convergence:
1810 #
1811 norm_u_new = util.Lsup(u_new)
1812 p_new=p+self.relaxation*g_new
1813 norm_p_new = util.sqrt(self.inner(p_new,p_new))
1814 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))
1815
1816 if self.iter>1:
1817 dg2=g_new-g
1818 df2=f_new-f
1819 norm_dg2=util.sqrt(self.inner(dg2,dg2))
1820 norm_df2=util.Lsup(df2)
1821 # print norm_g_new, norm_g, norm_dg, norm_p, tolerance
1822 tol_eq_g=tolerance*norm_dg2/(norm_g*abs(self.relaxation))*norm_p_new
1823 tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new
1824 if norm_g_new <= tol_eq_g and norm_f_new <= tol_eq_f:
1825 converged=True
1826 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
1827 self.trace("convergence after %s steps."%self.iter)
1828 return u,p

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26