/[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 2251 - (hide annotations)
Fri Feb 6 06:50:39 2009 UTC (11 years ago) by gross
File MIME type: text/x-python
File size: 64029 byte(s)
new compressible flow solver
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 2156 def GMRES(r, Aprod, x, bilinearform, atol=0, rtol=1.e-8, iter_max=100, iter_restart=20, verbose=False):
821     """
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 2245 x,stopped=_GMRESm(r2, Aprod, x, bilinearform, atol2, iter_max=iter_max-iter, iter_restart=m, verbose=verbose)
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 2245 def _GMRESm(r, Aprod, x, bilinearform, atol, iter_max=100, iter_restart=20, verbose=False):
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 2100 p=Aprod(v[iter])
914 artak 1519 v.append(p)
915    
916 caltinay 2169 v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
917 artak 1519
918 caltinay 2169 # Modified Gram-Schmidt
919     for j in range(iter+1):
920     h[j,iter]=bilinearform(v[j],v[iter+1])
921 gross 2100 v[iter+1]-=h[j,iter]*v[j]
922 caltinay 2169
923     h[iter+1,iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
924 gross 2100 v_norm2=h[iter+1,iter]
925 artak 1519
926     # Reorthogonalize if needed
927     if v_norm1 + 0.001*v_norm2 == v_norm1: #Brown/Hindmarsh condition (default)
928 caltinay 2169 for j in range(iter+1):
929 artak 1519 hr=bilinearform(v[j],v[iter+1])
930 caltinay 2169 h[j,iter]=h[j,iter]+hr
931 gross 1878 v[iter+1] -= hr*v[j]
932 artak 1519
933 caltinay 2169 v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
934 gross 2100 h[iter+1,iter]=v_norm2
935 artak 1519
936 caltinay 2169 # watch out for happy breakdown
937 gross 1878 if not v_norm2 == 0:
938 gross 2100 v[iter+1]=v[iter+1]/h[iter+1,iter]
939 artak 1519
940     # Form and store the information for the new Givens rotation
941 gross 2100 if iter > 0: h[:iter+1,iter]=__givapp(c[:iter],s[:iter],h[:iter+1,iter])
942     mu=math.sqrt(h[iter,iter]*h[iter,iter]+h[iter+1,iter]*h[iter+1,iter])
943 artak 1519
944     if mu!=0 :
945 gross 2100 c[iter]=h[iter,iter]/mu
946     s[iter]=-h[iter+1,iter]/mu
947     h[iter,iter]=c[iter]*h[iter,iter]-s[iter]*h[iter+1,iter]
948     h[iter+1,iter]=0.0
949 gross 1878 g[iter:iter+2]=__givapp(c[iter],s[iter],g[iter:iter+2])
950 artak 1519 # Update the residual norm
951 caltinay 2169
952 artak 1519 rho=abs(g[iter+1])
953 gross 2100 if verbose: print "GMRES: iteration step %s: residual %e"%(iter,rho)
954 artak 1519 iter+=1
955    
956 gross 1878 # At this point either iter > iter_max or rho < tol.
957 caltinay 2169 # It's time to compute x and leave.
958 gross 1878
959 gross 2100 if verbose: print "GMRES: iteration stopped after %s step."%iter
960 caltinay 2169 if iter > 0 :
961     y=numarray.zeros(iter,numarray.Float64)
962 gross 2100 y[iter-1] = g[iter-1] / h[iter-1,iter-1]
963 caltinay 2169 if iter > 1 :
964     i=iter-2
965 artak 1519 while i>=0 :
966 gross 2100 y[i] = ( g[i] - numarray.dot(h[i,i+1:iter], y[i+1:iter])) / h[i,i]
967 artak 1519 i=i-1
968     xhat=v[iter-1]*y[iter-1]
969     for i in range(iter-1):
970     xhat += v[i]*y[i]
971 caltinay 2169 else:
972 gross 2100 xhat=v[0] * 0
973 gross 1878
974 artak 1519 x += xhat
975 caltinay 2169 if iter<iter_restart-1:
976     stopped=True
977     else:
978 artak 1519 stopped=False
979    
980     return x,stopped
981    
982 gross 2156 def MINRES(r, Aprod, x, Msolve, bilinearform, atol=0, rtol=1.e-8, iter_max=100):
983     """
984 caltinay 2169 Solver for
985    
986 gross 2156 M{Ax=b}
987 caltinay 2169
988     with a symmetric and positive definite operator A (more details required!).
989     It uses the minimum residual method (MINRES) with preconditioner M
990     providing an approximation of A.
991    
992     The iteration is terminated if
993    
994 gross 2156 M{|r| <= atol+rtol*|r0|}
995 caltinay 2169
996 gross 2156 where M{r0} is the initial residual and M{|.|} is the energy norm. In fact
997    
998 caltinay 2169 M{|r| = sqrt( bilinearform(Msolve(r),r))}
999    
1000 gross 2156 For details on the preconditioned conjugate gradient method see the book:
1001    
1002 caltinay 2169 I{Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
1003     T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
1004     C. Romine, and H. van der Vorst}.
1005    
1006 gross 2156 @param r: initial residual M{r=b-Ax}. C{r} is altered.
1007     @type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1008 caltinay 2169 @param x: an initial guess for the solution
1009 gross 2156 @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1010     @param Aprod: returns the value Ax
1011 caltinay 2169 @type Aprod: function C{Aprod(x)} where C{x} is of the same object like
1012     argument C{x}. The returned object needs to be of the same
1013     type like argument C{r}.
1014 gross 2156 @param Msolve: solves Mx=r
1015 caltinay 2169 @type Msolve: function C{Msolve(r)} where C{r} is of the same type like
1016     argument C{r}. The returned object needs to be of the same
1017     type like argument C{x}.
1018 gross 2156 @param bilinearform: inner product C{<x,r>}
1019 caltinay 2169 @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same
1020     type like argument C{x} and C{r} is. The returned value
1021     is a C{float}.
1022 gross 2156 @param atol: absolute tolerance
1023     @type atol: non-negative C{float}
1024     @param rtol: relative tolerance
1025     @type rtol: non-negative C{float}
1026 caltinay 2169 @param iter_max: maximum number of iteration steps
1027 gross 2156 @type iter_max: C{int}
1028     @return: the solution approximation and the corresponding residual
1029 caltinay 2169 @rtype: C{tuple}
1030 gross 2156 @warning: C{r} and C{x} are altered.
1031     """
1032 artak 1481 #------------------------------------------------------------------
1033     # Set up y and v for the first Lanczos vector v1.
1034     # y = beta1 P' v1, where P = C**(-1).
1035     # v is really P' v1.
1036     #------------------------------------------------------------------
1037 gross 2156 r1 = r
1038     y = Msolve(r)
1039     beta1 = bilinearform(y,r)
1040 caltinay 2169
1041 artak 1481 if beta1< 0: raise NegativeNorm,"negative norm."
1042    
1043 gross 2156 # If r = 0 exactly, stop with x
1044     if beta1==0: return x
1045 artak 1481
1046 caltinay 2169 if beta1> 0: beta1 = math.sqrt(beta1)
1047 artak 1481
1048     #------------------------------------------------------------------
1049 artak 1484 # Initialize quantities.
1050 artak 1481 # ------------------------------------------------------------------
1051 artak 1482 iter = 0
1052     Anorm = 0
1053     ynorm = 0
1054 artak 1481 oldb = 0
1055     beta = beta1
1056     dbar = 0
1057     epsln = 0
1058     phibar = beta1
1059     rhs1 = beta1
1060     rhs2 = 0
1061     rnorm = phibar
1062     tnorm2 = 0
1063     ynorm2 = 0
1064     cs = -1
1065     sn = 0
1066 gross 2156 w = r*0.
1067     w2 = r*0.
1068 artak 1481 r2 = r1
1069     eps = 0.0001
1070    
1071     #---------------------------------------------------------------------
1072     # Main iteration loop.
1073     # --------------------------------------------------------------------
1074 gross 2100 while not rnorm<=atol+rtol*Anorm*ynorm: # checks ||r|| < (||A|| ||x||) * TOL
1075 artak 1481
1076     if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
1077     iter = iter + 1
1078    
1079     #-----------------------------------------------------------------
1080     # Obtain quantities for the next Lanczos vector vk+1, k = 1, 2,...
1081     # The general iteration is similar to the case k = 1 with v0 = 0:
1082     #
1083     # p1 = Operator * v1 - beta1 * v0,
1084     # alpha1 = v1'p1,
1085     # q2 = p2 - alpha1 * v1,
1086     # beta2^2 = q2'q2,
1087     # v2 = (1/beta2) q2.
1088     #
1089     # Again, y = betak P vk, where P = C**(-1).
1090     #-----------------------------------------------------------------
1091     s = 1/beta # Normalize previous vector (in y).
1092     v = s*y # v = vk if P = I
1093 caltinay 2169
1094 artak 1481 y = Aprod(v)
1095 caltinay 2169
1096 artak 1481 if iter >= 2:
1097     y = y - (beta/oldb)*r1
1098    
1099     alfa = bilinearform(v,y) # alphak
1100 caltinay 2169 y += (- alfa/beta)*r2
1101 artak 1481 r1 = r2
1102     r2 = y
1103     y = Msolve(r2)
1104     oldb = beta # oldb = betak
1105 artak 1550 beta = bilinearform(y,r2) # beta = betak+1^2
1106 artak 1481 if beta < 0: raise NegativeNorm,"negative norm."
1107    
1108     beta = math.sqrt( beta )
1109     tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta
1110 caltinay 2169
1111 artak 1481 if iter==1: # Initialize a few things.
1112     gmax = abs( alfa ) # alpha1
1113     gmin = gmax # alpha1
1114    
1115     # Apply previous rotation Qk-1 to get
1116     # [deltak epslnk+1] = [cs sn][dbark 0 ]
1117     # [gbar k dbar k+1] [sn -cs][alfak betak+1].
1118 caltinay 2169
1119 artak 1481 oldeps = epsln
1120     delta = cs * dbar + sn * alfa # delta1 = 0 deltak
1121     gbar = sn * dbar - cs * alfa # gbar 1 = alfa1 gbar k
1122     epsln = sn * beta # epsln2 = 0 epslnk+1
1123     dbar = - cs * beta # dbar 2 = beta2 dbar k+1
1124    
1125     # Compute the next plane rotation Qk
1126    
1127     gamma = math.sqrt(gbar*gbar+beta*beta) # gammak
1128 caltinay 2169 gamma = max(gamma,eps)
1129 artak 1481 cs = gbar / gamma # ck
1130     sn = beta / gamma # sk
1131     phi = cs * phibar # phik
1132     phibar = sn * phibar # phibark+1
1133    
1134     # Update x.
1135    
1136 caltinay 2169 denom = 1/gamma
1137     w1 = w2
1138     w2 = w
1139 artak 1481 w = (v - oldeps*w1 - delta*w2) * denom
1140 artak 1550 x += phi*w
1141 artak 1481
1142     # Go round again.
1143    
1144     gmax = max(gmax,gamma)
1145     gmin = min(gmin,gamma)
1146     z = rhs1 / gamma
1147     ynorm2 = z*z + ynorm2
1148     rhs1 = rhs2 - delta*z
1149     rhs2 = - epsln*z
1150    
1151     # Estimate various norms and test for convergence.
1152    
1153 caltinay 2169 Anorm = math.sqrt( tnorm2 )
1154     ynorm = math.sqrt( ynorm2 )
1155 artak 1481
1156     rnorm = phibar
1157    
1158     return x
1159 artak 1489
1160 gross 2156 def TFQMR(r, Aprod, x, bilinearform, atol=0, rtol=1.e-8, iter_max=100):
1161     """
1162 caltinay 2169 Solver for
1163 artak 1489
1164 gross 2156 M{Ax=b}
1165 artak 1489
1166 caltinay 2169 with a general operator A (more details required!).
1167     It uses the Transpose-Free Quasi-Minimal Residual method (TFQMR).
1168 artak 1489
1169 caltinay 2169 The iteration is terminated if
1170 artak 1489
1171 gross 2156 M{|r| <= atol+rtol*|r0|}
1172    
1173     where M{r0} is the initial residual and M{|.|} is the energy norm. In fact
1174    
1175 caltinay 2169 M{|r| = sqrt( bilinearform(r,r))}
1176 gross 2156
1177     @param r: initial residual M{r=b-Ax}. C{r} is altered.
1178     @type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1179 caltinay 2169 @param x: an initial guess for the solution
1180 gross 2156 @type x: same like C{r}
1181     @param Aprod: returns the value Ax
1182 caltinay 2169 @type Aprod: function C{Aprod(x)} where C{x} is of the same object like
1183     argument C{x}. The returned object needs to be of the same type
1184     like argument C{r}.
1185 gross 2156 @param bilinearform: inner product C{<x,r>}
1186 caltinay 2169 @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same
1187     type like argument C{x} and C{r}. The returned value is
1188     a C{float}.
1189 gross 2156 @param atol: absolute tolerance
1190     @type atol: non-negative C{float}
1191     @param rtol: relative tolerance
1192     @type rtol: non-negative C{float}
1193 caltinay 2169 @param iter_max: maximum number of iteration steps
1194 gross 2156 @type iter_max: C{int}
1195 caltinay 2169 @rtype: C{tuple}
1196 gross 2156 @warning: C{r} and C{x} are altered.
1197     """
1198 artak 1489 u1=0
1199     u2=0
1200     y1=0
1201     y2=0
1202    
1203     w = r
1204 caltinay 2169 y1 = r
1205     iter = 0
1206 artak 1489 d = 0
1207 gross 2156 v = Aprod(y1)
1208 artak 1489 u1 = v
1209 caltinay 2169
1210 artak 1489 theta = 0.0;
1211     eta = 0.0;
1212 gross 2156 rho=bilinearform(r,r)
1213     if rho < 0: raise NegativeNorm,"negative norm."
1214     tau = math.sqrt(rho)
1215     norm_r0=tau
1216     while tau>atol+rtol*norm_r0:
1217 artak 1489 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
1218    
1219     sigma = bilinearform(r,v)
1220 gross 2156 if sigma == 0.0: raise IterationBreakDown,'TFQMR breakdown, sigma=0'
1221 caltinay 2169
1222 artak 1489 alpha = rho / sigma
1223    
1224     for j in range(2):
1225     #
1226     # Compute y2 and u2 only if you have to
1227     #
1228     if ( j == 1 ):
1229     y2 = y1 - alpha * v
1230 gross 2156 u2 = Aprod(y2)
1231 caltinay 2169
1232 artak 1489 m = 2 * (iter+1) - 2 + (j+1)
1233 caltinay 2169 if j==0:
1234 artak 1489 w = w - alpha * u1
1235     d = y1 + ( theta * theta * eta / alpha ) * d
1236     if j==1:
1237     w = w - alpha * u2
1238     d = y2 + ( theta * theta * eta / alpha ) * d
1239    
1240     theta = math.sqrt(bilinearform(w,w))/ tau
1241     c = 1.0 / math.sqrt ( 1.0 + theta * theta )
1242     tau = tau * theta * c
1243     eta = c * c * alpha
1244     x = x + eta * d
1245     #
1246     # Try to terminate the iteration at each pass through the loop
1247     #
1248 gross 2156 if rho == 0.0: raise IterationBreakDown,'TFQMR breakdown, rho=0'
1249 artak 1489
1250     rhon = bilinearform(r,w)
1251     beta = rhon / rho;
1252     rho = rhon;
1253     y1 = w + beta * y2;
1254 gross 2156 u1 = Aprod(y1)
1255 artak 1489 v = u1 + beta * ( u2 + beta * v )
1256 caltinay 2169
1257 gross 2156 iter += 1
1258 artak 1489
1259     return x
1260    
1261    
1262 artak 1465 #############################################
1263    
1264 gross 1331 class ArithmeticTuple(object):
1265     """
1266 caltinay 2169 Tuple supporting inplace update x+=y and scaling x=a*y where C{x,y} is an
1267     ArithmeticTuple and C{a} is a float.
1268 gross 1331
1269 caltinay 2169 Example of usage::
1270 gross 1331
1271 caltinay 2169 from esys.escript import Data
1272     from numarray import array
1273     a=Data(...)
1274     b=array([1.,4.])
1275     x=ArithmeticTuple(a,b)
1276     y=5.*x
1277 gross 1331
1278     """
1279     def __init__(self,*args):
1280     """
1281 caltinay 2169 Initializes object with elements C{args}.
1282 gross 1331
1283 caltinay 2169 @param args: tuple of objects that support inplace add (x+=y) and
1284     scaling (x=a*y)
1285 gross 1331 """
1286     self.__items=list(args)
1287    
1288     def __len__(self):
1289     """
1290 caltinay 2169 Returns the number of items.
1291 gross 1331
1292     @return: number of items
1293     @rtype: C{int}
1294     """
1295     return len(self.__items)
1296    
1297     def __getitem__(self,index):
1298     """
1299 caltinay 2169 Returns item at specified position.
1300 gross 1331
1301 caltinay 2169 @param index: index of item to be returned
1302 gross 1331 @type index: C{int}
1303     @return: item with index C{index}
1304     """
1305     return self.__items.__getitem__(index)
1306    
1307     def __mul__(self,other):
1308     """
1309 caltinay 2169 Scales by C{other} from the right.
1310 gross 1331
1311     @param other: scaling factor
1312     @type other: C{float}
1313     @return: itemwise self*other
1314     @rtype: L{ArithmeticTuple}
1315     """
1316     out=[]
1317 caltinay 2169 try:
1318 gross 1896 l=len(other)
1319     if l!=len(self):
1320 caltinay 2169 raise ValueError,"length of arguments don't match."
1321 gross 1896 for i in range(l): out.append(self[i]*other[i])
1322     except TypeError:
1323 caltinay 2169 for i in range(len(self)): out.append(self[i]*other)
1324 gross 1331 return ArithmeticTuple(*tuple(out))
1325    
1326     def __rmul__(self,other):
1327     """
1328 caltinay 2169 Scales by C{other} from the left.
1329 gross 1331
1330     @param other: scaling factor
1331     @type other: C{float}
1332     @return: itemwise other*self
1333     @rtype: L{ArithmeticTuple}
1334     """
1335     out=[]
1336 caltinay 2169 try:
1337 gross 1896 l=len(other)
1338     if l!=len(self):
1339 caltinay 2169 raise ValueError,"length of arguments don't match."
1340 gross 1896 for i in range(l): out.append(other[i]*self[i])
1341     except TypeError:
1342 caltinay 2169 for i in range(len(self)): out.append(other*self[i])
1343 gross 1331 return ArithmeticTuple(*tuple(out))
1344    
1345 artak 1465 def __div__(self,other):
1346     """
1347 caltinay 2169 Scales by (1/C{other}) from the right.
1348 artak 1465
1349     @param other: scaling factor
1350     @type other: C{float}
1351     @return: itemwise self/other
1352     @rtype: L{ArithmeticTuple}
1353     """
1354 gross 1896 return self*(1/other)
1355 artak 1465
1356     def __rdiv__(self,other):
1357     """
1358 caltinay 2169 Scales by (1/C{other}) from the left.
1359 artak 1465
1360     @param other: scaling factor
1361     @type other: C{float}
1362     @return: itemwise other/self
1363     @rtype: L{ArithmeticTuple}
1364     """
1365     out=[]
1366 caltinay 2169 try:
1367 gross 1896 l=len(other)
1368     if l!=len(self):
1369 caltinay 2169 raise ValueError,"length of arguments don't match."
1370 gross 1896 for i in range(l): out.append(other[i]/self[i])
1371     except TypeError:
1372 caltinay 2169 for i in range(len(self)): out.append(other/self[i])
1373 artak 1465 return ArithmeticTuple(*tuple(out))
1374 caltinay 2169
1375 gross 1331 def __iadd__(self,other):
1376     """
1377 caltinay 2169 Inplace addition of C{other} to self.
1378 gross 1331
1379     @param other: increment
1380     @type other: C{ArithmeticTuple}
1381     """
1382     if len(self) != len(other):
1383 caltinay 2169 raise ValueError,"tuple lengths must match."
1384 gross 1331 for i in range(len(self)):
1385     self.__items[i]+=other[i]
1386     return self
1387    
1388 artak 1550 def __add__(self,other):
1389     """
1390 caltinay 2169 Adds C{other} to self.
1391 artak 1550
1392     @param other: increment
1393     @type other: C{ArithmeticTuple}
1394     """
1395 gross 1896 out=[]
1396 caltinay 2169 try:
1397 gross 1896 l=len(other)
1398     if l!=len(self):
1399 caltinay 2169 raise ValueError,"length of arguments don't match."
1400 gross 1896 for i in range(l): out.append(self[i]+other[i])
1401     except TypeError:
1402 caltinay 2169 for i in range(len(self)): out.append(self[i]+other)
1403 gross 1896 return ArithmeticTuple(*tuple(out))
1404 artak 1550
1405     def __sub__(self,other):
1406     """
1407 caltinay 2169 Subtracts C{other} from self.
1408 artak 1550
1409 caltinay 2169 @param other: decrement
1410 artak 1550 @type other: C{ArithmeticTuple}
1411     """
1412 gross 1896 out=[]
1413 caltinay 2169 try:
1414 gross 1896 l=len(other)
1415     if l!=len(self):
1416 caltinay 2169 raise ValueError,"length of arguments don't match."
1417 gross 1896 for i in range(l): out.append(self[i]-other[i])
1418     except TypeError:
1419 caltinay 2169 for i in range(len(self)): out.append(self[i]-other)
1420 gross 1896 return ArithmeticTuple(*tuple(out))
1421 caltinay 2169
1422 artak 1557 def __isub__(self,other):
1423     """
1424 caltinay 2169 Inplace subtraction of C{other} from self.
1425 artak 1550
1426 caltinay 2169 @param other: decrement
1427 artak 1557 @type other: C{ArithmeticTuple}
1428     """
1429     if len(self) != len(other):
1430     raise ValueError,"tuple length must match."
1431     for i in range(len(self)):
1432     self.__items[i]-=other[i]
1433     return self
1434    
1435 artak 1550 def __neg__(self):
1436     """
1437 caltinay 2169 Negates values.
1438 artak 1550 """
1439 gross 1896 out=[]
1440 artak 1550 for i in range(len(self)):
1441 gross 1896 out.append(-self[i])
1442     return ArithmeticTuple(*tuple(out))
1443 artak 1550
1444    
1445 gross 1414 class HomogeneousSaddlePointProblem(object):
1446     """
1447 caltinay 2169 This class provides a framework for solving linear homogeneous saddle
1448     point problems of the form::
1449 gross 1414
1450 caltinay 2169 M{Av+B^*p=f}
1451     M{Bv =0}
1452 gross 1414
1453 caltinay 2169 for the unknowns M{v} and M{p} and given operators M{A} and M{B} and
1454     given right hand side M{f}. M{B^*} is the adjoint operator of M{B}.
1455 gross 1414 """
1456     def __init__(self,**kwargs):
1457     self.setTolerance()
1458 gross 2100 self.setAbsoluteTolerance()
1459 gross 2156 self.setSubProblemTolerance()
1460 gross 1414
1461 gross 2100 #=============================================================
1462 gross 1414 def initialize(self):
1463     """
1464 caltinay 2169 Initializes the problem (overwrite).
1465 gross 1414 """
1466     pass
1467 gross 2100
1468 gross 2251 def inner_pBv(self,p,v):
1469 gross 1414 """
1470 caltinay 2169 Returns inner product of element p and Bv (overwrite).
1471    
1472 gross 2251 @param p: a pressure increment
1473     @param v: a residual
1474     @return: inner product of element p and Bv
1475     @rtype: C{float}
1476     @note: used if PCG is applied.
1477 gross 1414 """
1478 gross 2100 raise NotImplementedError,"no inner product for p implemented."
1479 gross 1414
1480 gross 2100 def inner_p(self,p0,p1):
1481 gross 1414 """
1482 gross 2251 Returns inner product of p0 and p1 (overwrite).
1483 caltinay 2169
1484 gross 2251 @param p0: a pressure
1485     @param p1: a pressure
1486     @return: inner product of p0 and p1
1487     @rtype: C{float}
1488 gross 1414 """
1489 gross 2100 raise NotImplementedError,"no inner product for p implemented."
1490 gross 2251
1491     def norm_v(self,v):
1492     """
1493     Returns the norm of v (overwrite).
1494 gross 2100
1495 gross 2251 @param v: a velovity
1496     @return: norm of v
1497     @rtype: non-negative C{float}
1498 gross 2100 """
1499 gross 2251 raise NotImplementedError,"no norm of v implemented."
1500 caltinay 2169
1501 gross 2251
1502     def getV(self, p, v0):
1503 gross 2100 """
1504 gross 2251 return the value for v for a given p (overwrite)
1505 gross 1414
1506 gross 2251 @param p: a pressure
1507     @param v0: a initial guess for the value v to return.
1508     @return: v given as M{v= A^{-1} (f-B^*p)}
1509 gross 1414 """
1510 gross 2251 raise NotImplementedError,"no v calculation implemented."
1511    
1512     def norm_Bv(self,v):
1513     """
1514     Returns Bv (overwrite).
1515    
1516     @rtype: equal to the type of p
1517     @note: boundary conditions on p should be zero!
1518     """
1519     raise NotImplementedError, "no operator B implemented."
1520    
1521     def solve_AinvBt(self,p):
1522     """
1523     Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}
1524 caltinay 2169 (overwrite).
1525 gross 1414
1526 gross 2251 @param p: a pressure increment
1527     @return: the solution of M{Av=B^*p}
1528 gross 1414 @note: boundary conditions on v should be zero!
1529     """
1530 gross 2100 raise NotImplementedError,"no operator A implemented."
1531 gross 1414
1532 gross 2251 def solve_precB(self,v):
1533 gross 1414 """
1534 caltinay 2169 Provides a preconditioner for M{BA^{-1}B^*} with accuracy
1535     L{self.getSubProblemTolerance()} (overwrite).
1536 gross 1414
1537     @rtype: equal to the type of p
1538 gross 2100 @note: boundary conditions on p should be zero!
1539 gross 1414 """
1540 gross 2100 raise NotImplementedError,"no preconditioner for Schur complement implemented."
1541     #=============================================================
1542     def __Aprod_PCG(self,p):
1543 gross 2251 return self.solve_AinvBt(p)
1544 gross 1414
1545 gross 2251 def __inner_PCG(self,p,v):
1546     return self.inner_pBv(p,v)
1547 gross 1414
1548 gross 2251 def __Msolve_PCG(self,v):
1549     return self.solve_precB(v)
1550 gross 2100 #=============================================================
1551 gross 2251 def __Aprod_GMRES(self,p):
1552     return self.solve_precB(self.solve_AinvBt(p))
1553     def __inner_GMRES(self,p0,p1):
1554     return self.inner_p(p0,p1)
1555 gross 2100 #=============================================================
1556 gross 2251 def norm_p(self,p):
1557     """
1558     calculates the norm of C{p}
1559    
1560     @param p: a pressure
1561     @return: the norm of C{p} using the inner product for pressure
1562     @rtype: C{float}
1563     """
1564     f=self.inner_p(p,p)
1565     if f<0: raise ValueError,"negative pressure norm."
1566     return math.sqrt(f)
1567    
1568 artak 1550
1569 gross 2251 def solve(self,v,p,max_iter=20, verbose=False, show_details=False, usePCG=True, iter_restart=20, max_correction_steps=10):
1570 gross 2100 """
1571 caltinay 2169 Solves the saddle point problem using initial guesses v and p.
1572 gross 1414
1573 gross 2100 @param v: initial guess for velocity
1574     @param p: initial guess for pressure
1575     @type v: L{Data}
1576     @type p: L{Data}
1577 gross 2251 @param usePCG: indicates the usage of the PCG rather than GMRES scheme.
1578 caltinay 2169 @param max_iter: maximum number of iteration steps per correction
1579     attempt
1580     @param verbose: if True, shows information on the progress of the
1581     saddlepoint problem solver.
1582     @param show_details: if True, shows details of the sub problem solver
1583     @param iter_restart: restart the iteration after C{iter_restart} steps
1584     (only used if useUzaw=False)
1585 gross 2251 @type usePCG: C{bool}
1586 gross 2100 @type max_iter: C{int}
1587     @type verbose: C{bool}
1588     @type show_details: C{bool}
1589     @type iter_restart: C{int}
1590     @rtype: C{tuple} of L{Data} objects
1591     """
1592     self.verbose=verbose
1593     self.show_details=show_details and self.verbose
1594     rtol=self.getTolerance()
1595     atol=self.getAbsoluteTolerance()
1596     correction_step=0
1597     converged=False
1598 gross 2251 while not converged:
1599     # calculate velocity for current pressure:
1600     v=self.getV(p,v)
1601     #
1602     norm_v=self.norm_v(v)
1603     norm_Bv=self.norm_Bv(v)
1604     ATOL=norm_v*rtol+atol
1605     if self.verbose: print "saddle point solver: norm v= %e, norm_Bv= %e, tolerance = %e."%(norm_v, norm_Bv,ATOL)
1606     if not ATOL>0: raise ValueError,"overall absolute tolerance needs to be positive."
1607     if norm_Bv <= ATOL:
1608     converged=True
1609     else:
1610     correction_step+=1
1611     if correction_step>max_correction_steps:
1612     raise CorrectionFailed,"Given up after %d correction steps."%correction_step
1613     dp=self.solve_precB(v)
1614     if usePCG:
1615     norm2=self.inner_pBv(dp,v)
1616     if norm2<0: raise ValueError,"negative PCG norm."
1617     norm2=math.sqrt(norm2)
1618     else:
1619     norm2=self.norm_p(dp)
1620     ATOL_ITER=ATOL/norm_Bv*norm2
1621     if self.verbose: print "saddle point solver: tolerance for solver: %e"%ATOL_ITER
1622     if usePCG:
1623     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)
1624     else:
1625     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)
1626 gross 2100 if self.verbose: print "saddle point solver: tolerance reached."
1627 gross 2251 return v,p
1628 artak 1465
1629 caltinay 2169 #========================================================================
1630 gross 2100 def setTolerance(self,tolerance=1.e-4):
1631     """
1632 caltinay 2169 Sets the relative tolerance for (v,p).
1633 gross 2100
1634 caltinay 2169 @param tolerance: tolerance to be used
1635 gross 2100 @type tolerance: non-negative C{float}
1636     """
1637     if tolerance<0:
1638     raise ValueError,"tolerance must be positive."
1639     self.__rtol=tolerance
1640 gross 2156 self.setSubProblemTolerance()
1641    
1642 gross 1414 def getTolerance(self):
1643 gross 2100 """
1644 caltinay 2169 Returns the relative tolerance.
1645 gross 1414
1646 caltinay 2169 @return: relative tolerance
1647 gross 2100 @rtype: C{float}
1648     """
1649     return self.__rtol
1650 caltinay 2169
1651 gross 2100 def setAbsoluteTolerance(self,tolerance=0.):
1652     """
1653 caltinay 2169 Sets the absolute tolerance.
1654 gross 1414
1655 caltinay 2169 @param tolerance: tolerance to be used
1656 gross 2100 @type tolerance: non-negative C{float}
1657     """
1658     if tolerance<0:
1659     raise ValueError,"tolerance must be non-negative."
1660     self.__atol=tolerance
1661 caltinay 2169
1662 gross 2100 def getAbsoluteTolerance(self):
1663     """
1664 caltinay 2169 Returns the absolute tolerance.
1665 gross 1414
1666 caltinay 2169 @return: absolute tolerance
1667 gross 2100 @rtype: C{float}
1668     """
1669     return self.__atol
1670 gross 1469
1671 gross 2156 def setSubProblemTolerance(self,rtol=None):
1672 gross 2100 """
1673 caltinay 2169 Sets the relative tolerance to solve the subproblem(s).
1674 artak 1489
1675 gross 2100 @param rtol: relative tolerence
1676     @type rtol: positive C{float}
1677     """
1678 caltinay 2169 if rtol == None:
1679 gross 2156 rtol=max(200.*util.EPSILON,self.getTolerance()**2)
1680 gross 2100 if rtol<=0:
1681     raise ValueError,"tolerance must be positive."
1682 gross 2156 self.__sub_tol=rtol
1683 caltinay 2169
1684 gross 2100 def getSubProblemTolerance(self):
1685     """
1686 caltinay 2169 Returns the subproblem reduction factor.
1687 artak 1557
1688 gross 2100 @return: subproblem reduction factor
1689 caltinay 2169 @rtype: C{float}
1690 gross 2100 """
1691     return self.__sub_tol
1692 artak 1550
1693 gross 1878 def MaskFromBoundaryTag(domain,*tags):
1694     """
1695 caltinay 2169 Creates a mask on the Solution(domain) function space where the value is
1696     one for samples that touch the boundary tagged by tags.
1697 gross 1878
1698 caltinay 2169 Usage: m=MaskFromBoundaryTag(domain, "left", "right")
1699 gross 1878
1700 caltinay 2169 @param domain: domain to be used
1701 gross 1878 @type domain: L{escript.Domain}
1702 caltinay 2169 @param tags: boundary tags
1703 gross 1878 @type tags: C{str}
1704 caltinay 2169 @return: a mask which marks samples that are touching the boundary tagged
1705     by any of the given tags
1706 gross 1878 @rtype: L{escript.Data} of rank 0
1707     """
1708     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1709 gross 1956 d=escript.Scalar(0.,escript.FunctionOnBoundary(domain))
1710 gross 1878 for t in tags: d.setTaggedValue(t,1.)
1711     pde.setValue(y=d)
1712     return util.whereNonZero(pde.getRightHandSide())
1713 caltinay 2169
1714     #==============================================================================
1715 gross 867 class SaddlePointProblem(object):
1716     """
1717 caltinay 2169 This implements a solver for a saddle point problem
1718 gross 867
1719 gross 877 M{f(u,p)=0}
1720     M{g(u)=0}
1721 gross 867
1722     for u and p. The problem is solved with an inexact Uszawa scheme for p:
1723    
1724 ksteube 990 M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})}
1725 gross 877 M{Q_g (p^{k+1}-p^{k}) = g(u^{k+1})}
1726 gross 867
1727 caltinay 2169 where Q_f is an approximation of the Jacobian A_f of f with respect to u and
1728     Q_f is an approximation of A_g A_f^{-1} A_g with A_g is the Jacobian of g
1729     with respect to p. As a the construction of a 'proper' Q_g can be difficult,
1730     non-linear conjugate gradient method is applied to solve for p, so Q_g plays
1731 gross 867 in fact the role of a preconditioner.
1732     """
1733     def __init__(self,verbose=False,*args):
1734     """
1735 caltinay 2169 Initializes the problem.
1736 gross 867
1737 caltinay 2169 @param verbose: if True, some information is printed in the process
1738 gross 867 @type verbose: C{bool}
1739 caltinay 2169 @note: this method may be overwritten by a particular saddle point
1740     problem
1741 gross 867 """
1742 gross 1659 print "SaddlePointProblem should not be used anymore!"
1743 gross 1107 if not isinstance(verbose,bool):
1744     raise TypeError("verbose needs to be of type bool.")
1745 gross 1106 self.__verbose=verbose
1746 gross 877 self.relaxation=1.
1747 gross 1878 DeprecationWarning("SaddlePointProblem should not be used anymore.")
1748 gross 867
1749     def trace(self,text):
1750     """
1751 caltinay 2169 Prints C{text} if verbose has been set.
1752 gross 867
1753 ksteube 990 @param text: a text message
1754 gross 867 @type text: C{str}
1755     """
1756 ksteube 1567 if self.__verbose: print "%s: %s"%(str(self),text)
1757 gross 867
1758 gross 873 def solve_f(self,u,p,tol=1.e-8):
1759 gross 867 """
1760 caltinay 2169 Solves
1761 gross 867
1762 caltinay 2169 A_f du = f(u,p)
1763 gross 867
1764 caltinay 2169 with tolerance C{tol} and returns du. A_f is Jacobian of f with respect
1765     to u.
1766 gross 867
1767     @param u: current approximation of u
1768     @type u: L{escript.Data}
1769     @param p: current approximation of p
1770     @type p: L{escript.Data}
1771 gross 873 @param tol: tolerance expected for du
1772 gross 867 @type tol: C{float}
1773     @return: increment du
1774     @rtype: L{escript.Data}
1775 caltinay 2169 @note: this method has to be overwritten by a particular saddle point
1776     problem
1777 gross 867 """
1778     pass
1779    
1780 gross 873 def solve_g(self,u,tol=1.e-8):
1781 gross 867 """
1782 caltinay 2169 Solves
1783 gross 867
1784 caltinay 2169 Q_g dp = g(u)
1785 gross 867
1786 caltinay 2169 where Q_g is a preconditioner for A_g A_f^{-1} A_g with A_g is the
1787     Jacobian of g with respect to p.
1788 gross 867
1789     @param u: current approximation of u
1790     @type u: L{escript.Data}
1791 gross 873 @param tol: tolerance expected for dp
1792     @type tol: C{float}
1793 gross 867 @return: increment dp
1794     @rtype: L{escript.Data}
1795 caltinay 2169 @note: this method has to be overwritten by a particular saddle point
1796     problem
1797 gross 867 """
1798     pass
1799    
1800     def inner(self,p0,p1):
1801     """
1802 caltinay 2169 Inner product of p0 and p1 approximating p. Typically this returns
1803     C{integrate(p0*p1)}.
1804 gross 867 @return: inner product of p0 and p1
1805     @rtype: C{float}
1806     """
1807     pass
1808    
1809 gross 877 subiter_max=3
1810     def solve(self,u0,p0,tolerance=1.e-6,tolerance_u=None,iter_max=100,accepted_reduction=0.995,relaxation=None):
1811     """
1812 caltinay 2169 Runs the solver.
1813 gross 873
1814 gross 877 @param u0: initial guess for C{u}
1815     @type u0: L{esys.escript.Data}
1816     @param p0: initial guess for C{p}
1817     @type p0: L{esys.escript.Data}
1818     @param tolerance: tolerance for relative error in C{u} and C{p}
1819     @type tolerance: positive C{float}
1820 caltinay 2169 @param tolerance_u: tolerance for relative error in C{u} if different
1821     from C{tolerance}
1822 gross 877 @type tolerance_u: positive C{float}
1823 caltinay 2169 @param iter_max: maximum number of iteration steps
1824 gross 877 @type iter_max: C{int}
1825 caltinay 2169 @param accepted_reduction: if the norm g cannot be reduced by
1826     C{accepted_reduction} backtracking to adapt
1827     the relaxation factor. If
1828     C{accepted_reduction=None} no backtracking
1829     is used.
1830 gross 877 @type accepted_reduction: positive C{float} or C{None}
1831 caltinay 2169 @param relaxation: initial relaxation factor. If C{relaxation==None},
1832     the last relaxation factor is used.
1833 gross 877 @type relaxation: C{float} or C{None}
1834     """
1835     tol=1.e-2
1836     if tolerance_u==None: tolerance_u=tolerance
1837     if not relaxation==None: self.relaxation=relaxation
1838     if accepted_reduction ==None:
1839     angle_limit=0.
1840     elif accepted_reduction>=1.:
1841     angle_limit=0.
1842     else:
1843     angle_limit=util.sqrt(1-accepted_reduction**2)
1844     self.iter=0
1845     u=u0
1846     p=p0
1847     #
1848     # initialize things:
1849     #
1850     converged=False
1851     #
1852     # start loop:
1853     #
1854     # initial search direction is g
1855     #
1856     while not converged :
1857     if self.iter>iter_max:
1858     raise ArithmeticError("no convergence after %s steps."%self.iter)
1859     f_new=self.solve_f(u,p,tol)
1860     norm_f_new = util.Lsup(f_new)
1861     u_new=u-f_new
1862     g_new=self.solve_g(u_new,tol)
1863     self.iter+=1
1864     norm_g_new = util.sqrt(self.inner(g_new,g_new))
1865     if norm_f_new==0. and norm_g_new==0.: return u, p
1866     if self.iter>1 and not accepted_reduction==None:
1867     #
1868     # did we manage to reduce the norm of G? I
1869     # if not we start a backtracking procedure
1870     #
1871     # print "new/old norm = ",norm_g_new, norm_g, norm_g_new/norm_g
1872     if norm_g_new > accepted_reduction * norm_g:
1873     sub_iter=0
1874     s=self.relaxation
1875     d=g
1876     g_last=g
1877     self.trace(" start substepping: f = %s, g = %s, relaxation = %s."%(norm_f_new, norm_g_new, s))
1878     while sub_iter < self.subiter_max and norm_g_new > accepted_reduction * norm_g:
1879     dg= g_new-g_last
1880     norm_dg=abs(util.sqrt(self.inner(dg,dg))/self.relaxation)
1881     rad=self.inner(g_new,dg)/self.relaxation
1882     # print " ",sub_iter,": rad, norm_dg:",abs(rad), norm_dg*norm_g_new * angle_limit
1883     # print " ",sub_iter,": rad, norm_dg:",rad, norm_dg, norm_g_new, norm_g
1884     if abs(rad) < norm_dg*norm_g_new * angle_limit:
1885     if sub_iter>0: self.trace(" no further improvements expected from backtracking.")
1886     break
1887     r=self.relaxation
1888     self.relaxation= - rad/norm_dg**2
1889     s+=self.relaxation
1890     #####
1891     # a=g_new+self.relaxation*dg/r
1892     # print "predicted new norm = ",util.sqrt(self.inner(a,a)),util.sqrt(self.inner(g_new,g_new)), self.relaxation
1893     #####
1894     g_last=g_new
1895     p+=self.relaxation*d
1896     f_new=self.solve_f(u,p,tol)
1897     u_new=u-f_new
1898     g_new=self.solve_g(u_new,tol)
1899     self.iter+=1
1900     norm_f_new = util.Lsup(f_new)
1901     norm_g_new = util.sqrt(self.inner(g_new,g_new))
1902     # print " ",sub_iter," new g norm",norm_g_new
1903     self.trace(" %s th sub-step: f = %s, g = %s, relaxation = %s."%(sub_iter, norm_f_new, norm_g_new, s))
1904     #
1905     # can we expect reduction of g?
1906     #
1907     # u_last=u_new
1908     sub_iter+=1
1909     self.relaxation=s
1910     #
1911     # check for convergence:
1912     #
1913     norm_u_new = util.Lsup(u_new)
1914     p_new=p+self.relaxation*g_new
1915     norm_p_new = util.sqrt(self.inner(p_new,p_new))
1916 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))
1917 gross 873
1918 gross 877 if self.iter>1:
1919     dg2=g_new-g
1920     df2=f_new-f
1921     norm_dg2=util.sqrt(self.inner(dg2,dg2))
1922     norm_df2=util.Lsup(df2)
1923     # print norm_g_new, norm_g, norm_dg, norm_p, tolerance
1924     tol_eq_g=tolerance*norm_dg2/(norm_g*abs(self.relaxation))*norm_p_new
1925     tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new
1926     if norm_g_new <= tol_eq_g and norm_f_new <= tol_eq_f:
1927     converged=True
1928     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
1929     self.trace("convergence after %s steps."%self.iter)
1930     return u,p
1931 caltinay 2169

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26