/[escript]/trunk/escript/py_src/pdetools.py
ViewVC logotype

Annotation of /trunk/escript/py_src/pdetools.py

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26