/[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 2455 - (hide annotations)
Wed Jun 3 03:29:07 2009 UTC (10 years, 8 months ago) by jfenwick
File MIME type: text/x-python
File size: 64950 byte(s)
Merging changes from numpy branch.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26