/[escript]/branches/arrexp_trunk2098/escript/py_src/pdetools.py
ViewVC logotype

Contents of /branches/arrexp_trunk2098/escript/py_src/pdetools.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2145 - (show annotations)
Wed Dec 10 00:48:53 2008 UTC (10 years, 6 months ago) by jfenwick
File MIME type: text/x-python
File size: 59564 byte(s)
Branch commit.
This version of the python files passes unit tests compiled with 
usenumarray=no and usenumarray=yes
That is, it does not use any functions returning 
boost::python::numeric::array


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26