/[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 1878 - (show annotations)
Tue Oct 14 03:39:13 2008 UTC (10 years, 11 months ago) by gross
File MIME type: text/x-python
File size: 59595 byte(s)
new version of JacobiFree Newton GMRES + test added.
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 other=1.*other
1261 if isinstance(other,float):
1262 for i in range(len(self)):
1263 out.append(self[i]*other)
1264 else:
1265 for i in range(len(self)):
1266 out.append(self[i]*other[i])
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 other=1.*other
1280 if isinstance(other,float):
1281 for i in range(len(self)):
1282 out.append(other*self[i])
1283 else:
1284 for i in range(len(self)):
1285 out.append(other[i]*self[i])
1286 return ArithmeticTuple(*tuple(out))
1287
1288 #########################
1289 # Added by Artak
1290 #########################
1291 def __div__(self,other):
1292 """
1293 dividing from the right
1294
1295 @param other: scaling factor
1296 @type other: C{float}
1297 @return: itemwise self/other
1298 @rtype: L{ArithmeticTuple}
1299 """
1300 out=[]
1301 other=1.*other
1302 if isinstance(other,float):
1303 for i in range(len(self)):
1304 out.append(self[i]*(1./other))
1305 else:
1306 for i in range(len(self)):
1307 out.append(self[i]*(1./other[i]))
1308 return ArithmeticTuple(*tuple(out))
1309
1310 def __rdiv__(self,other):
1311 """
1312 dividing from the left
1313
1314 @param other: scaling factor
1315 @type other: C{float}
1316 @return: itemwise other/self
1317 @rtype: L{ArithmeticTuple}
1318 """
1319 out=[]
1320 other=1.*other
1321 if isinstance(other,float):
1322 for i in range(len(self)):
1323 out.append(other*(1./self[i]))
1324 else:
1325 for i in range(len(self)):
1326 out.append(other[i]*(1./self[i]))
1327 return ArithmeticTuple(*tuple(out))
1328
1329 ##########################################33
1330
1331 def __iadd__(self,other):
1332 """
1333 in-place add of other to self
1334
1335 @param other: increment
1336 @type other: C{ArithmeticTuple}
1337 """
1338 if len(self) != len(other):
1339 raise ValueError,"tuple length must match."
1340 for i in range(len(self)):
1341 self.__items[i]+=other[i]
1342 return self
1343
1344 def __add__(self,other):
1345 """
1346 add other to self
1347
1348 @param other: increment
1349 @type other: C{ArithmeticTuple}
1350 """
1351 # if not isinstance(other,float):
1352 if len(self) != len(other):
1353 raise ValueError,"tuple length must match."
1354 for i in range(len(self)):
1355 self.__items[i]+=other[i]
1356 # else:
1357 # for i in range(len(self)):
1358 # self.__items[i]+=other
1359
1360 return self
1361
1362 def __sub__(self,other):
1363 """
1364 subtract other from self
1365
1366 @param other: increment
1367 @type other: C{ArithmeticTuple}
1368 """
1369 if len(self) != len(other):
1370 raise ValueError,"tuple length must match."
1371 for i in range(len(self)):
1372 self.__items[i]-=other[i]
1373 return self
1374
1375 def __isub__(self,other):
1376 """
1377 subtract other from self
1378
1379 @param other: increment
1380 @type other: C{ArithmeticTuple}
1381 """
1382 if len(self) != len(other):
1383 raise ValueError,"tuple length must match."
1384 for i in range(len(self)):
1385 self.__items[i]-=other[i]
1386 return self
1387
1388 def __neg__(self):
1389 """
1390 negate
1391
1392 """
1393 for i in range(len(self)):
1394 self.__items[i]=-self.__items[i]
1395 return self
1396
1397
1398 class HomogeneousSaddlePointProblem(object):
1399 """
1400 This provides a framwork for solving linear homogeneous saddle point problem of the form
1401
1402 Av+B^*p=f
1403 Bv =0
1404
1405 for the unknowns v and p and given operators A and B and given right hand side f.
1406 B^* is the adjoint operator of B is the given inner product.
1407
1408 """
1409 def __init__(self,**kwargs):
1410 self.setTolerance()
1411 self.setToleranceReductionFactor()
1412
1413 def initialize(self):
1414 """
1415 initialize the problem (overwrite)
1416 """
1417 pass
1418 def B(self,v):
1419 """
1420 returns Bv (overwrite)
1421 @rtype: equal to the type of p
1422
1423 @note: boundary conditions on p should be zero!
1424 """
1425 pass
1426
1427 def inner(self,p0,p1):
1428 """
1429 returns inner product of two element p0 and p1 (overwrite)
1430
1431 @type p0: equal to the type of p
1432 @type p1: equal to the type of p
1433 @rtype: C{float}
1434
1435 @rtype: equal to the type of p
1436 """
1437 pass
1438
1439 def solve_A(self,u,p):
1440 """
1441 solves Av=f-Au-B^*p with accuracy self.getReducedTolerance() (overwrite)
1442
1443 @rtype: equal to the type of v
1444 @note: boundary conditions on v should be zero!
1445 """
1446 pass
1447
1448 def solve_prec(self,p):
1449 """
1450 provides a preconditioner for BA^{-1}B^* with accuracy self.getReducedTolerance() (overwrite)
1451
1452 @rtype: equal to the type of p
1453 """
1454 pass
1455
1456 def stoppingcriterium(self,Bv,v,p):
1457 """
1458 returns a True if iteration is terminated. (overwrite)
1459
1460 @rtype: C{bool}
1461 """
1462 pass
1463
1464 def __inner(self,p,r):
1465 return self.inner(p,r[1])
1466
1467 def __inner_p(self,p1,p2):
1468 return self.inner(p1,p2)
1469
1470 def __inner_a(self,a1,a2):
1471 return self.inner_a(a1,a2)
1472
1473 def __inner_a1(self,a1,a2):
1474 return self.inner(a1[1],a2[1])
1475
1476 def __stoppingcriterium(self,norm_r,r,p):
1477 return self.stoppingcriterium(r[1],r[0],p)
1478
1479 def __stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None):
1480 return self.stoppingcriterium2(norm_r,norm_b,solver,TOL)
1481
1482 def setTolerance(self,tolerance=1.e-8):
1483 self.__tol=tolerance
1484 def getTolerance(self):
1485 return self.__tol
1486 def setToleranceReductionFactor(self,reduction=0.01):
1487 self.__reduction=reduction
1488 def getSubProblemTolerance(self):
1489 return self.__reduction*self.getTolerance()
1490
1491 def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='PCG',iter_restart=20):
1492 """
1493 solves the saddle point problem using initial guesses v and p.
1494
1495 @param max_iter: maximum number of iteration steps.
1496 """
1497 self.verbose=verbose
1498 self.show_details=show_details and self.verbose
1499
1500 # assume p is known: then v=A^-1(f-B^*p)
1501 # which leads to BA^-1B^*p = BA^-1f
1502
1503 # Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask)
1504 self.__z=v+self.solve_A(v,p*0)
1505 Bz=self.B(self.__z)
1506 #
1507 # solve BA^-1B^*p = Bz
1508 #
1509 #
1510 #
1511 self.iter=0
1512 if solver=='GMRES':
1513 if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter
1514 p=GMRES(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.,iter_restart=iter_restart)
1515 # solve Au=f-B^*p
1516 # A(u-v)=f-B^*p-Av
1517 # u=v+(u-v)
1518 u=v+self.solve_A(v,p)
1519
1520 if solver=='TFQMR':
1521 if self.verbose: print "enter TFQMR method (iter_max=%s)"%max_iter
1522 p=TFQMR(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)
1523 # solve Au=f-B^*p
1524 # A(u-v)=f-B^*p-Av
1525 # u=v+(u-v)
1526 u=v+self.solve_A(v,p)
1527
1528 if solver=='MINRES':
1529 if self.verbose: print "enter MINRES method (iter_max=%s)"%max_iter
1530 p=MINRES(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)
1531 # solve Au=f-B^*p
1532 # A(u-v)=f-B^*p-Av
1533 # u=v+(u-v)
1534 u=v+self.solve_A(v,p)
1535
1536 if solver=='GMRESC':
1537 if self.verbose: print "enter GMRES coupled method (iter_max=%s)"%max_iter
1538 p0=self.solve_prec1(Bz)
1539 #alfa=(1/self.vol)*util.integrate(util.interpolate(p,escript.Function(self.domain)))
1540 #p-=alfa
1541 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)
1542 #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())
1543
1544 # solve Au=f-B^*p
1545 # A(u-v)=f-B^*p-Av
1546 # u=v+(u-v)
1547 p=x[1]
1548 u=v+self.solve_A(v,p)
1549 #p=x[1]
1550 #u=x[0]
1551
1552 if solver=='PCG':
1553 # note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv
1554 #
1555 # with Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask)
1556 # A(v-z)= f -Az - B^*p (v-z=0 on fixed_u_mask)
1557 if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter
1558 p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p)
1559 u=r[0]
1560 # print "DIFF=",util.integrate(p)
1561
1562 # print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)
1563
1564 return u,p
1565
1566 def __Msolve(self,r):
1567 return self.solve_prec(r[1])
1568
1569 def __Msolve2(self,r):
1570 return self.solve_prec(r*1)
1571
1572 def __Mempty(self,r):
1573 return r
1574
1575
1576 def __Aprod(self,p):
1577 # return BA^-1B*p
1578 #solve Av =B^*p as Av =f-Az-B^*(-p)
1579 v=self.solve_A(self.__z,-p)
1580 return ArithmeticTuple(v, self.B(v))
1581
1582 def __Anext(self,x):
1583 # return next v,p
1584 #solve Av =-B^*p as Av =f-Az-B^*p
1585
1586 pc=x[1]
1587 v=self.solve_A(self.__z,-pc)
1588 p=self.solve_prec1(self.B(v))
1589
1590 return ArithmeticTuple(v,p)
1591
1592
1593 def __Aprod2(self,p):
1594 # return BA^-1B*p
1595 #solve Av =B^*p as Av =f-Az-B^*(-p)
1596 v=self.solve_A(self.__z,-p)
1597 return self.B(v)
1598
1599 def __Aprod_Newton(self,p):
1600 # return BA^-1B*p - Bz
1601 #solve Av =-B^*p as Av =f-Az-B^*p
1602 v=self.solve_A(self.__z,-p)
1603 return self.B(v-self.__z)
1604
1605 def __Aprod_Newton2(self,x):
1606 # return BA^-1B*p - Bz
1607 #solve Av =-B^*p as Av =f-Az-B^*p
1608 pc=x[1]
1609 v=self.solve_A(self.__z,-pc)
1610 p=self.solve_prec1(self.B(v-self.__z))
1611 return ArithmeticTuple(v,p)
1612
1613
1614 def MaskFromBoundaryTag(domain,*tags):
1615 """
1616 create a mask on the Solution(domain) function space which one for samples
1617 that touch the boundary tagged by tags.
1618
1619 usage: m=MaskFromBoundaryTag(domain,"left", "right")
1620
1621 @param domain: a given domain
1622 @type domain: L{escript.Domain}
1623 @param tags: boundray tags
1624 @type tags: C{str}
1625 @return: a mask which marks samples that are touching the boundary tagged by any of the given tags.
1626 @rtype: L{escript.Data} of rank 0
1627 """
1628 pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1629 d=escript.Scalar(0.,escript.FunctionOnBoundary(function_space.getDomain()))
1630 for t in tags: d.setTaggedValue(t,1.)
1631 pde.setValue(y=d)
1632 return util.whereNonZero(pde.getRightHandSide())
1633 #============================================================================================================================================
1634 class SaddlePointProblem(object):
1635 """
1636 This implements a solver for a saddlepoint problem
1637
1638 M{f(u,p)=0}
1639 M{g(u)=0}
1640
1641 for u and p. The problem is solved with an inexact Uszawa scheme for p:
1642
1643 M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})}
1644 M{Q_g (p^{k+1}-p^{k}) = g(u^{k+1})}
1645
1646 where Q_f is an approximation of the Jacobiean A_f of f with respect to u and Q_f is an approximation of
1647 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'
1648 Q_g can be difficult, non-linear conjugate gradient method is applied to solve for p, so Q_g plays
1649 in fact the role of a preconditioner.
1650 """
1651 def __init__(self,verbose=False,*args):
1652 """
1653 initializes the problem
1654
1655 @param verbose: switches on the printing out some information
1656 @type verbose: C{bool}
1657 @note: this method may be overwritten by a particular saddle point problem
1658 """
1659 print "SaddlePointProblem should not be used anymore!"
1660 if not isinstance(verbose,bool):
1661 raise TypeError("verbose needs to be of type bool.")
1662 self.__verbose=verbose
1663 self.relaxation=1.
1664 DeprecationWarning("SaddlePointProblem should not be used anymore.")
1665
1666 def trace(self,text):
1667 """
1668 prints text if verbose has been set
1669
1670 @param text: a text message
1671 @type text: C{str}
1672 """
1673 if self.__verbose: print "%s: %s"%(str(self),text)
1674
1675 def solve_f(self,u,p,tol=1.e-8):
1676 """
1677 solves
1678
1679 A_f du = f(u,p)
1680
1681 with tolerance C{tol} and return du. A_f is Jacobiean of f with respect to u.
1682
1683 @param u: current approximation of u
1684 @type u: L{escript.Data}
1685 @param p: current approximation of p
1686 @type p: L{escript.Data}
1687 @param tol: tolerance expected for du
1688 @type tol: C{float}
1689 @return: increment du
1690 @rtype: L{escript.Data}
1691 @note: this method has to be overwritten by a particular saddle point problem
1692 """
1693 pass
1694
1695 def solve_g(self,u,tol=1.e-8):
1696 """
1697 solves
1698
1699 Q_g dp = g(u)
1700
1701 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.
1702
1703 @param u: current approximation of u
1704 @type u: L{escript.Data}
1705 @param tol: tolerance expected for dp
1706 @type tol: C{float}
1707 @return: increment dp
1708 @rtype: L{escript.Data}
1709 @note: this method has to be overwritten by a particular saddle point problem
1710 """
1711 pass
1712
1713 def inner(self,p0,p1):
1714 """
1715 inner product of p0 and p1 approximating p. Typically this returns integrate(p0*p1)
1716 @return: inner product of p0 and p1
1717 @rtype: C{float}
1718 """
1719 pass
1720
1721 subiter_max=3
1722 def solve(self,u0,p0,tolerance=1.e-6,tolerance_u=None,iter_max=100,accepted_reduction=0.995,relaxation=None):
1723 """
1724 runs the solver
1725
1726 @param u0: initial guess for C{u}
1727 @type u0: L{esys.escript.Data}
1728 @param p0: initial guess for C{p}
1729 @type p0: L{esys.escript.Data}
1730 @param tolerance: tolerance for relative error in C{u} and C{p}
1731 @type tolerance: positive C{float}
1732 @param tolerance_u: tolerance for relative error in C{u} if different from C{tolerance}
1733 @type tolerance_u: positive C{float}
1734 @param iter_max: maximum number of iteration steps.
1735 @type iter_max: C{int}
1736 @param accepted_reduction: if the norm g cannot be reduced by C{accepted_reduction} backtracking to adapt the
1737 relaxation factor. If C{accepted_reduction=None} no backtracking is used.
1738 @type accepted_reduction: positive C{float} or C{None}
1739 @param relaxation: initial relaxation factor. If C{relaxation==None}, the last relaxation factor is used.
1740 @type relaxation: C{float} or C{None}
1741 """
1742 tol=1.e-2
1743 if tolerance_u==None: tolerance_u=tolerance
1744 if not relaxation==None: self.relaxation=relaxation
1745 if accepted_reduction ==None:
1746 angle_limit=0.
1747 elif accepted_reduction>=1.:
1748 angle_limit=0.
1749 else:
1750 angle_limit=util.sqrt(1-accepted_reduction**2)
1751 self.iter=0
1752 u=u0
1753 p=p0
1754 #
1755 # initialize things:
1756 #
1757 converged=False
1758 #
1759 # start loop:
1760 #
1761 # initial search direction is g
1762 #
1763 while not converged :
1764 if self.iter>iter_max:
1765 raise ArithmeticError("no convergence after %s steps."%self.iter)
1766 f_new=self.solve_f(u,p,tol)
1767 norm_f_new = util.Lsup(f_new)
1768 u_new=u-f_new
1769 g_new=self.solve_g(u_new,tol)
1770 self.iter+=1
1771 norm_g_new = util.sqrt(self.inner(g_new,g_new))
1772 if norm_f_new==0. and norm_g_new==0.: return u, p
1773 if self.iter>1 and not accepted_reduction==None:
1774 #
1775 # did we manage to reduce the norm of G? I
1776 # if not we start a backtracking procedure
1777 #
1778 # print "new/old norm = ",norm_g_new, norm_g, norm_g_new/norm_g
1779 if norm_g_new > accepted_reduction * norm_g:
1780 sub_iter=0
1781 s=self.relaxation
1782 d=g
1783 g_last=g
1784 self.trace(" start substepping: f = %s, g = %s, relaxation = %s."%(norm_f_new, norm_g_new, s))
1785 while sub_iter < self.subiter_max and norm_g_new > accepted_reduction * norm_g:
1786 dg= g_new-g_last
1787 norm_dg=abs(util.sqrt(self.inner(dg,dg))/self.relaxation)
1788 rad=self.inner(g_new,dg)/self.relaxation
1789 # print " ",sub_iter,": rad, norm_dg:",abs(rad), norm_dg*norm_g_new * angle_limit
1790 # print " ",sub_iter,": rad, norm_dg:",rad, norm_dg, norm_g_new, norm_g
1791 if abs(rad) < norm_dg*norm_g_new * angle_limit:
1792 if sub_iter>0: self.trace(" no further improvements expected from backtracking.")
1793 break
1794 r=self.relaxation
1795 self.relaxation= - rad/norm_dg**2
1796 s+=self.relaxation
1797 #####
1798 # a=g_new+self.relaxation*dg/r
1799 # print "predicted new norm = ",util.sqrt(self.inner(a,a)),util.sqrt(self.inner(g_new,g_new)), self.relaxation
1800 #####
1801 g_last=g_new
1802 p+=self.relaxation*d
1803 f_new=self.solve_f(u,p,tol)
1804 u_new=u-f_new
1805 g_new=self.solve_g(u_new,tol)
1806 self.iter+=1
1807 norm_f_new = util.Lsup(f_new)
1808 norm_g_new = util.sqrt(self.inner(g_new,g_new))
1809 # print " ",sub_iter," new g norm",norm_g_new
1810 self.trace(" %s th sub-step: f = %s, g = %s, relaxation = %s."%(sub_iter, norm_f_new, norm_g_new, s))
1811 #
1812 # can we expect reduction of g?
1813 #
1814 # u_last=u_new
1815 sub_iter+=1
1816 self.relaxation=s
1817 #
1818 # check for convergence:
1819 #
1820 norm_u_new = util.Lsup(u_new)
1821 p_new=p+self.relaxation*g_new
1822 norm_p_new = util.sqrt(self.inner(p_new,p_new))
1823 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))
1824
1825 if self.iter>1:
1826 dg2=g_new-g
1827 df2=f_new-f
1828 norm_dg2=util.sqrt(self.inner(dg2,dg2))
1829 norm_df2=util.Lsup(df2)
1830 # print norm_g_new, norm_g, norm_dg, norm_p, tolerance
1831 tol_eq_g=tolerance*norm_dg2/(norm_g*abs(self.relaxation))*norm_p_new
1832 tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new
1833 if norm_g_new <= tol_eq_g and norm_f_new <= tol_eq_f:
1834 converged=True
1835 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
1836 self.trace("convergence after %s steps."%self.iter)
1837 return u,p

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26